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Introduction 


In  this  work,  a  novel  imaging  technique  is  explored  that  uses  non-harmful  application  of  near  infrared  light  to  determine  the 
pathophysiological  properties  of  breast  tissue.  Using  this  technique,  known  as  near  infrared  (NIR)  tomography,  an  optical  fiber  placed 
on  the  surface  of  the  region  of  interest,  the  breast,  delivers  an  input  signal  while  other  optical  fibers  placed  at  different  locations  on  the 
same  surface  detect  the  out-coming  photons,  which  have  propagated  through  the  volume  under  investigation.  The  intensity  and  path- 
length  distributions  of  the  exiting  light  provide  information  about  the  optical  properties  of  the  transilluminated  tissue  using  a  model- 
based  interpretation  where  photon  propagation  is  simulated  by  the  diffusion  theory.  Through  iterative  solution  to  match  the  theory  to 
the  real  measurements,  images  of  internal  absorption  and  scattering  coefficient  distribution  can  be  reconstructed.  Tissue  optical 
properties  have  shown  to  be  a  function  of  their  structure  and  more  importantly,  their  physiological  state.  It  is  these  qualities  that 
provide  an  alternative  method  for  the  detection  and  characterization  of  tumor  within  breast  tissue.  An  important  aspect  for  the  success 
and  accuracy  of  this  method  is  the  use  of  adequate  modeling  of  light  within  breast  tissue.  In  this  work,  we  have  used  a  Finite  Element 
Model  of  the  Diffusion  Approximation,  to  calculate  and  predict  light  transport.  Using  such  technique,  we  are  able  to  determine,  to 
some  accuracy,  the  distribution  of  internal  optical  properties  of  the  tissue  under  investigation.  A  major  limitation  in  the  development 
of  near  infrared  imaging  so  far  has  been  the  lack  of  accurate  and  fast  model-based  reconstruction  algorithm  for  this  modality,  which  is 
a  direct  aim  of  the  project  being  funded. 

Since  the  award  of  this  project,  the  aims  have  not  changed.  In  brief  these  are: 

(1)  Validate,  improve  and  accelerate  a  three-dimensional  finite  element  model  of  light  propagation  within  the  breast  tissue  as 
well  as  the  image  reconstruction  algorithm  for  simultaneously  calculating  the  internal  absorption  and  scattering  coefficients. 

(2)  Investigate  the  benefits  and  limits  of  using  a-priori  information  from  dual  modality  images,  for  example,  MRI  and  Near 
Infrared  data. 

(3)  Explore  the  benefits  and  limits  of  using  Fluorescent  contrast-agents,  which  will  target  specific  molecular  markers  of  the 
tumor,  giving  rise  to  a-priori  information  regarding  tumor  location  as  well  as  physiological  function  of  the  tissue. 

(4)  Test  the  developed  algorithm  on  patient  data  and  compare  results  with  available  mammograms  and  biopsy  results  to  evaluate 
algorithm  accuracy 

During  the  second  year  of  the  funding  period,  several  advances  and  improvements  to  the  modeling  and  image  reconstruction 
algorithms  have  been  implemented.  The  three  dimensional  image  reconstruction  algorithms  have  been  further  modified  to  improve 
computational  speed  and  overhead  allowing  a  more  realistic  time  frame  for  image  reconstruction.  Additional  improvements  in  the 
modeling  capabilities  have  been  incorporated,  namely  in  terms  of  allowing  for  internal  tissue  Refractive  Index  (RI)  distribution  as  well 
as  tissue  deformation,  both  of  which  are  discussed  below.  The  use  of  a-priori  information  within  the  image  reconstruction  algorithm 
have  been  investigates,  namely  through  the  use  of  structural  (MRI)  information  and  spectral  (Wavelength  dependant)  information. 
Finally  some  image  analysis  has  been  used  to  date,  utilizing  the  developed  code  to  seek  the  pathophysiological  properties  of  breast 
tissue  and  cancer.  Each  of  the  area  of  research  developed  during  this  second  of  funding  will  be  addressed  below  together  with  the 
main  findings  and  their  implications  of  optical  tomography  of  breast  tissue  in-vivo,  specifically  for  the  detection  and  characterization 
of  breast  cancer. 


Effect  of  Refractive  Index  variation 

Previous  studies  [1]  reported  results  in  modeling  the  effect  of  Refractive  Index 
(RI)  on  the  forward  model.  In  the  results  obtained  during  this  funding  period,  the 
initiative  was  advanced  to  investigating  the  effect  of  RI  variation  on  NIR  image 
reconstruction,  by  assuming  either  correct  knowledge  of  the  RI  of  each  tissue  or 
applying  a  homogenous  value  throughout  the  model.  The  study  showed  that 
providing  the  RI  of  the  glandular  tissue  is  not  far  from  the  value  of  adipose  tissue, 
it  has  little  effect  on  the  qualitative  and  quantitative  accuracy  of  the  results,  as 
demonstrated  in  Figure  1.  Assuming  the  RI  of  the  whole  breast  is  similar,  for 
example,  the  RI  is  1.455  for  adipose  and  1.4  for  glandular  tissue,  reconstructed 
images  of  absorption  and  scatter  can  be  obtained  which  ignore  the  effect  of  RI 
with  modest  degradation  in  the  recovery  of  information  about  an  abnormality. 
However,  it  is  also  important  to  note  that  this  analysis  was  completed  under  the 
assumption  that  the  abnormality  (i.e.  tumor)  has  the  same  RI  as  its  background 
tissue  and  is  most  typically  located  in  the  fibroglandular  tissue.  Further  studies  are 
required  to  establish  the  RI  variation  for  different  types  of  tumors  and  to 
investigate  how  such  variations  might  alter  NIR  tomography.  The  absolute  values 
of  RI  for  breast  tissue  types  are  still  a  subject  under  investigation.  Estimates  are 
difficult  to  obtain  accurately  in-vivo,  but  some  data  has  been  reported  for  various 
tissue  types  [2,  3].  Although  adipose  tissue  (fatty  layer)  has  been  measured  to  be 
1.455,  no  results  exist  for  fibro-glandular  tissue.  Nonetheless  the  RI  of  glandular 
tissue  is  believed  to  be  lower  than  that  of  adipose  (e.g.  values  of  1.4  have  been 
assumed  in  [2,  4]).  The  rationale  for  a  lower  value  is  sound  in  that  estimates  of  the 
composition  of  fibroglandular  tissue  place  its  water  (typically  over  60%)  and 
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Figure  1.  Images  reconstructed  using  a  model  with 
refractive  index  variation.  Top  row  shows  the  exact 
optical  properties  where  adipose  and  glandular  tissue 
have  different  optical  properties  and  different  refractive 
index  (RI)  of  1.455  and  1.4,  respectively.  Middle  row 
shows  reconstructed  images  assuming  the  correct  RI 
distribution.  Bottom  row  shows  reconstructed  images 
assuming  a  homogenous  RI  of  1 .455. 


blood  content  (over  1%)  to  be  high,  indicating  that  its  bulk  refractive  index  is  likely  close  to  that  of  water  [5,  6]. 

Breast  deformation  modeling 

In  order  to  obtain 
good  NIR  data 
measurements  it 
is  essential  to 
have  good 

contact  between 
the  optical  fibers 
and  the  breast 
which  in-turn 
results  in  the 
deformation  of 
the  breast,  Figure 
2,  due  to  the  soft 

plasticity  of  the 
tissue.  A  tissue 
deformation 
model  of  the 
female  breast  has 
been  developed 
that  will  account 
for  the  altered 
shape  of  the 
breast  during 
clinical  NIR 
measurements, 

Figure  3.  Using  a 
deformed  model 
of  a  breast,  NIR 

data  was  used  to  reconstruct  images  of  tissue  absorption  and  reduced  scatter  using  assumptions  about  the  imaging  domain,  Figure  4 
[7], 

Assuming  a  non-deformed  breast  shape  for  image  reconstruction  has  shown  to  lead  to  poor  quality  images  since  the  geometry  of  the 
breast  is  greatly  altered,  whereas  using  the  correct  deformed  geometry  produces  the  best  images.  The  goal  of  this  work  was  to 
incorporate  this  new  model  of  breast  deformation  with  more  accurate  information  regarding  the  mechanical  properties  of  the  breast  to 
improve  the  NIR  image  reconstruction.  The  mechanical  property  information  is  readily  available  from  other  imaging  modalities,  and 
the  synthesis  of  this  information  may  provide  fundamentally  new  information  about  breast  physiologic  response  to  pressure,  and/or 
breast  pathology  response  to  pressure.  An  accurate  model  of  breast  deformation  should  in  principle  allow  us  to  create  patient  specific 
models  and  meshes,  which  would  in-turn  provide  more  clinically  useful  data. 

Spectral  and  spatial  a-priori  information 

The  use  of  spatial  a  priori  information  was  investigated  during  the  first  year  of  funding  under  this  program.  Namely,  realistic  breast 
phantoms  were  used  to  investigate  the  use  of  spatial  information  to  constrain  and  penalize  the  reconstruction  algorithm  and  have  been 
shown  to  improve  both  the  qualitative  and  quantitative  accuracy  of  the  reconstructed  images  [8-10]. 

The  concept  underpinning  parametric  or  spectrally  constrained  imaging  is  the  generation  of  a 
derived  response  from  multi- spectral  data.  Typically,  images  are  formed  serially,  wavelength- 
by-wavelength,  and  do  not  impose  connectivity  between  the  image  estimate  at  one 
wavelength  with  that  of  at  another.  By  using  a  spectral  model  composed  of  certain 
wavelength-independent  parameters,  as  in  Figure  5,  it  becomes  possible  to  combine  multi- 
spectral  data  in  a  new  way.  Concepts  of  multi-spectral  and  spatially  constrained  image 
reconstruction  have  recently  been  extended  to  NIR  which  has  reached  the  point  of  application 
to  experimental  and  clinical  data  and  is  highlighted  here  as  an  illustration  in  Figure  6  (row  3). 

Specifically,  multiple  wavelengths  NIR  data  have  been  used  simultaneously  to  reconstruct 
absolute  values  of  chromophore  concentration  (deoxy  and  oxygenated  blood  concentration 
and  water)  as  well  as  scattering  properties  (scattering  power  (b)  and  amplitude  (A)).  Thus, 
instead  of  estimating  the  absorption  and  reduced  scatter  properties  at  each  wavelength, 
independently  and  then  using  these  values  to  derive  specific  chromophore  concentrations  and 
scattering  parameters,  the  aim  is  to  encode  directly  the  linear  map  between  wavelength-dependent  optical  properties  and  chromophore 
concentrations  through  their  known  spectral  signatures.  This  leads  to  a  single  reconstruction  to  estimate  all  of  the  final  images 
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Figure  5.  Absorption  coefficients  of 
various  chromophores  as  a  function  of  NIR 
wavelengths. 


Figure  2.  Photo  of  the  breast  being 
compressed  by  the  optical  fibers  to  ensure 
good  contact  for  data  collection. 


Figure  3.  Volume  mesh  of  (a)  the  normal  suspended  breast  and  (b)  the  deformed  mesh 
after  the  application  of  the  optical  fiber  array. 


(b) 


Reduced  Scatter 


Reduced  Scatter 


Figure  4.  (a)  2D  coronal  slices  through  the  deformed  breast  mesh,  showing  the  position  of  the  anomalies.  The  most  right  hand  slice  is 
near  the  chest  while  the  most  left  hand  slice  is  near  the  nipple,  (b)  Reconstructed  images  assuming  no  deformation,  and  (c)  reconstructed 
images  using  a  deformed  mesh  using  information  about  fiber  contact. 
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simultaneously  and  has  shown  to  offer  the  advantage  of  reducing  image  sensitivity  to  noise  at  specific  (individual)  wavelengths, 
especially  in  the  water  and  oxygen  saturation  distributions  [11,  12]. 

A  combined  NIR-MRI  imaging  system  has  been  used  [9,  1 0]  in  a 
case  study  to  estimate  the  properties  of  healthy  breast  tissue, 

Figure  6.  The  images  obtained  for  the  NIR  parameters  using  an 
unconstrained  reconstruction  shows  noisy  images  with  boundary 
artifacts.  The  spatial  priors  act  on  these  images,  making  them 
smoother,  but  preserving  the  trend  for  the  quantification.  For 
example,  the  scatter  power  shows  a  reduction  in  the  glandular 
tissue  (row  2)  similar  to  the  values  obtained  in  the  no-priors 
reconstruction  (row  1).  However,  previous  studies  have  indicated 
that,  glandular  tissue  having  higher  number  density  of  scatterers, 
may  actually  have  higher  values  for  scatter  power,  than  fatty 
tissue.  Hence,  the  results  from  spatially  constrained 
reconstruction,  while  appearing  smoother,  may  be  misleading. 

The  scatter  power  image  obtained  by  application  of  the  spectrally 
constrained  method  is  more  quantitatively  acceptable  and  the 
spatial  priors  acting  on  this  spectral  method  gives  the  most 
intuitively  acceptable  image  for  this  parameter,  showing  the 
layered  structure  of  the  breast  tissue.  We  observed  elevated 
[HbT]  (25:13  jliM),  water  (91:49  %),  and  scattering  power 
(1.0:0. 5)  in  glandular  tissue  relative  to  adipose  tissue  using  the 
combined  priors,  which  matches  the  higher  vascularization  in  the 
glandular  tissue. 

The  results  so  far  show  that,  while  anatomical  information 
improves  the  image  quality  resulting  in  reduced  artifacts,  it  may 
not  significantly  improve  quantification.  The  spectral  prior 
obtained  using  the  intrinsic  behavior  of  tissue  chromophores  and 
scattering  in  the  near  infrared  wavelength  regime  plays  a  more 
important  role  in  addressing  this  problem,  and  finally,  a  synergy  between  these  two  priors  yields  the  most  accurate  characterization  of 
breast  tissue  properties  currently  possible. 

Tissue  characterization  using  Spectral  Imaging 

As  an  illustration,  the  direct  spectral  reconstruction  method  has  also  been  applied 
to  clinical  tomography  data  obtained  from  measurements  on  a  66  year  old  subject 
with  a  5-10  mm  spiculated  mass  evident  on  mammography,  which  was  later 
diagnosed  as  Invasive  Ductal  Carcinoma  [12].  The  subject  underwent  the  NIR 
exam  in  accordance  with  the  Dartmouth  protocol,  where  informed  consent  was 
obtained  prior  to  the  exam.  The  position  of  the  anomaly  was  provided  through 
mammography  and  determined  to  be  at  1 1:30  o’clock  in  the  cranocaudal  view  and 
located  about  1cm  from  surface.  Reconstructed  images  are  shown  in  Figure  7.  A 
localized  increase  in  [HbT]  was  observed  and  the  contrast  available  was  1.7: 1.0  in 
tumor  versus  background.  The  [St02]  image  showed  a  decrease  at  the  location  of 
the  tumor  with  a  contrast  of  0.7:1.  The  corresponding  St02  image  with  the 
conventional  technique  was  noisier  (spatially)  and  a  similar  decrease  was  not 
observed.  The  water  fraction  image  showed  heterogeneity  in  the  range  0.25  to 
0.5%,  which  is  considerably  smaller  than  the  range  of  0  to  1%  observed  with  the 
old  method.  Tighter  data  ranges  are  similarly  observed  in  the  scatter  images,  and 
artifacts  evident  in  the  separate  wavelength  technique  have  been  suppressed  in  the 
spectral  approach. 

With  this  new  technique,  higher  hemoglobin  content  at  the  site  of  tumor  is 
observed  in  the  clinical  case,  as  well  as  lower  oxygen  saturation.  The  water 
fraction  and  scattering  images  show  the  most  significant  improvement  in  all  cases  by  reduction  of  the  image  artifacts.  Further  data 
analysis  over  a  larger  patient  data  is  needed  and  is  a  subject  of  active  research  in  the  final  year  of  the  funded  project. 
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Figure  7.  HbT  [mM],  St02  [%],  water  fraction,  ‘a’  and 
‘b’  images  are  shown  from  measurements  on  a  cancer 
subject  with  a  5-10  mm  Invasive  Ductal  Carcinoma  at 
the  1 1 :30  clock  face  position. 
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Figure  6.  Breast  tissue  property  images  for  a  healthy  female  volunteer 
estimated  using  four  different  reconstruction  methods.  First  (top),  only  the 
outer  boundary  of  the  imaging  domain,  and  the  location  of  the  optical  fibre 
measurement  sites  are  specified.  Second,  a  spatially  constrained  algorithm  was 
used.  The  MRI  of  this  patient  defined  the  spatial  constraints,  which  relate  to 
the  internal  distribution  of  adipose  and  glandular  tissues.  Third,  spectral 
constraints  were  applied  and  chromophore  concentrations  and  scattering 
parameters  were  reconstructed  directly.  Fourth  (bottom),  both  spatial  and 
spectral  constraints  were  combined.  Reconstructed  images  are  for  total  blood 
(HbT),  Oxygen  saturation  (St02),  Water  (H20),  Scattering  amplitude  (A)  and 
Scattering  power  (b). 
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Key  research  accomplishments 


1.  Assuming  that  the  effect  of  the  refractive  index  (RI)  variation  in  breast  tissue  is  small,  such  that  if  ignored,  the  errors 
obtained  within  reconstructed  images  are  minimal.  This  in  reality  is  acceptable,  given  that  the  variation  in  RI  of  tissue  within 
the  breast  is  small,  ranging  from  1.4  to  1.455. 

2.  The  effect  of  breast  deformation  tissue  to  optical  fiber  contact  has  been  investigated.  It  is  shown  that  if  the  deformation  is 
large,  adequate  priori  knowledge  about  the  shape  of  the  breast  is  essential  for  useful  clinical  images.  Such  information  can 
either  be  obtained  using  mechanical  models  (as  demonstrated)  or  using  a-priori  structural  images  from  for  example  MRI. 

3.  Spectral  and  spatial  priori  constrains  have  been  developed  and  implemented  to  show  improved  accuracy  in  clinical 
information  about  reconstructed  images.  It  is  shown  that  incorporation  of  both  such  information  significantly  improves  both 
qualitative  and  quantitative  accuracy  of  images. 

4.  To  date,  spectral  constraints  have  been  used  to  analyze  clinical  data.  It  has  been  demonstrated  that  addition  of  such 
constraints  can  significantly  improve  the  clinical  accuracy  of  physiological  data  obtained  from  clinical  cases. 

Reportable  outcomes 

The  first  2  years  of  the  funding  of  this  project  has  led  to  a  number  of  peer  reviews  publications  as  well  as  oral  presentation  at 
international  conferences.  Links  and  collaborations  have  been  established  with  Washington  University  at  St.  Louis  who  are  now 
routinely  using  the  developed  code.  Further  collaboration  has  been  established  with  University  of  Pennsylvania  through  the 
development  of  alternative  image  reconstruction  algorithms  allowing  the  use  of  large  datasets.  Reportable  outcomes  include: 

1 .  Oral  presentation  at  the  Optical  Society  of  America  (OSA)  international  meeting  in  2004 

2.  Oral  presentation  at  the  International  Society  for  Optical  Engineering  international  meeting  in  2005 

3.  Applied  for  funding  from  NCI  through  R01  mechanism  (PI  Hamid  Dehghani,  “Breast  Deformation  Modeling  for  Near 
Infrared  Tomography”  3  year  project) 

4.  Applied  for  funding  from  DoD  through  the  IDEA  award  mechanism  (PI  Hamid  Dehghani,  “Neoadjuvant  Therapy  Treatment 
monitoring  using  CT  Guided  Optical  Tomography”  3  year  project) 

Conclusions 


The  project  is  continuing  to  reconstruct  NIR  images  of  clinical  data  to  build  the  database  of  3D  images.  Steps  have  also  been  taken  as 
part  of  this  project  to  achieve  other  aims  specified  by  this  project.  A  new  and  novel  image  reconstruction  algorithm  has  been 
developed  and  implemented  which  allows  the  direct  reconstruction  of  chromophore  concentration  and  scattering  functions  of  breast 
tissue.  This  not  only  allows  a  computationally  efficient  implementations,  but  also  since  all  the  data  from  all  wavelength  are  used 
coupled  together,  a  great  reduction  in  image  noise  and  artifact  is  seen.  Finally,  as  a  direct  result  of  the  2nd  aim,  the  project  is  beginning 
to  develop,  understand  and  implement  a  useful  algorithm  for  the  incorporation  of  a-priori  information,  not  only  in  phantom  studies, 
but  also  in  clinical  setups. 

References 


1.  Dehghani,  H.,  Brooksby,  B.,  Vishwanath,  K.,  Pogue,  B.  W.,  Paulsen  K.  D.,  The  effects  of  internal  refractive  index  variation 
in  near  infrared  optical  tomography:  A  finite  element  modeling  approach.  Phys.  Med.  Biol.,  2003.  48:  p.  2713-2727. 

2.  Bolin,  F.P.,  Preuss,  L.  E.,  Taylor,  R.  C.,  Ference,  R.  J,  Refractive  index  of  some  mammalian  tissue  using  a  fiber  optic 
cladding  method.  Applied  Optics,  1989.  28(12):  p.  2297-2303. 

3.  Brooksby,  B.,  Dehghani,  H.,  Vishwanath,  K.,  Pogue,  B.  W.,  Paulsen,  K.  D.,  Internal  refractive  index  changes  affect  light 
transport  in  tissue,  proceedings  of  spie,  2003.  4955:  p.  296-304. 

4.  Tromberg,  B.J.,  Coquoz,  O.,  Fishkin,  J.  B.,  Pham,  T.,  Anderson,  E.  R.,  Butler,  J.,  Cahn,  M.,  Gross,  J.  D.,  Venugopalan,  V., 
Pham,  D.,  Non-invasive  measurements  of  breast  tissue  optical  properties  using  frequency-domain  photon  migration.  Phil. 
Trans.  R.  Soc.  Lond.  B,  1997.  352:  p.  661-668. 

5.  Srinivasan,  S.,  Pogue,  B.  W.,  Jiang,  S.,  Dehghani,  H.,  Kogel,  C.,  Soho,  S.,  Gibson,  J.  J.,  Tosteson,  T.  D.,  Poplack,  S.  P., 
Paulsen,  K.  D.,  Interpreting  Hemoglobin  and  Water  Concentration,  Oxygen  Saturation  and  Scattering  Measured  In  Vivo  by 
Near-Infrared  Breast  Tomography.  PNAS,  2003.  100(21):  p.  12349-12354. 

6.  Pogue,  B.W.,  Jiang,  S.,  Dehghani,  H.,  Kogel,  C.,  Soho,  S.,  Srinivasan,  S.,  Song,  X.,  Poplack,  S.  P.,  Paulsen,  K.  D., 
Characterization  of  Hemoglobin,  Water  and  NIR  Scattering  in  Breast  Tissue:  Analysis  of  inter-subject  variability  and 
menstrual  cycle  changes  relative  to  lesions.  JBO,  2004.  9(3):  p.  541-552. 

7.  Dehghani,  H.,  Doyley,  M.  M.,  Pogue,  B.  W.,  Jiang,  S.,  Geng,  J.,  and  Paulsen,  K.  D.,  Breast  deformation  modeling  for  image 
reconstruction  in  near  infrared  optical  tomography,  phys.  med.  biol,  2004.  49(7):  p.  1131-1146. 

8.  Xu,  H.,  Springett,  R.,  Dehghani,  H.,  Pogue,  B.W.,  Paulsen,  K.D.,  and  Dunn,  J.F.,  MRI  coupled  Broadband  Near  Infrared 
Tomography  System  for  Small  Animal  Brain  Studies.  Appl.  Opt.,  2005.  44(1 1):  p.  2177-2188. 

9.  Brooksby,  B.,  Jiang,  S.,  Dehghani,  H.,  Pogue,  B.  W.,  Paulsen,  K.  D.,  Weaver,  J.,  Kogel,  C.,  and  Poplack,  S.  P.,  Combining 
near  infrared  tomography  and  magnetic  resonance  imaging  to  study  in  vivo  breast  tissue:  implementation  of  a  Laplacian- 
type  regularization  to  incorporate  magentic  reasonance  structure.  JBO,  2005.  10(5). 

7 


10.  Brooksby,  B.,  Jiang,  S.,  Kogel,  C.,  Doyley,  M.,  Dehghani,  H.,  Weaver,  J.  B.,  Poplack,  S.  P.,  Pogue,  B.  W.,  Paulsen,  K.  D., 
Magnetic  Resonance- Guided  Near-Infrared  Tomography  of  the  Breast.  Rev.  Sci.  Instrum.,  2004.  75(12). 

11.  Brooksby,  B.,  Srinivasan,  S.,  Jiang,  S.,  Dehghani,  H.,  Pogue,  B.W.,  Paulsen,  K.D.,  Weaver,  J.,  Kogel,  C.,  and  Poplack,  S.P., 
Spectral-prior  information  improves  Near-Infrared  diffuse  tomography  more  than  spatial-prior.  Optics  letters,  2005.  30(15): 
p.  1968-1970. 

12.  Srinivasan,  S.,  Pogue,  B.  W.,  Jiang,  S.,  Dehghani,  H.,  and  Paulsen,  K.  D.,  Spectrally  Constrained  Chromophore  and 
Scattering  NIR  Tomography  Provides  Quantitative  and  Robust  Reconstruction.  Appl.  Opt.,  2005.  44(10):  p.  1858-1869. 


Publications  List 

S.  Srinivasan,  B.  W.  Pogue,  S.  Jiang,  H.  Dehghani,  C.  Kogel,  S.  Soho,  J.  J.  Gibson,  T.  D.  Tosteson,  S.  P.  Poplack,  K.  D.  Paulsen, 
"Interpreting  Hemoglobin  and  Water  Concentration,  Oxygen  Saturation  and  Scattering  Measured  In  Vivo  by  Near-Infrared  Breast 
Tomography"  PNAS  vol  100  no  21  pp  12349-12354  2003 

X.  Song,  B.  W.  Pogue,  s.  Jiang,  M.  M.  Doyley,  H.  Dehghani,  T.  D.  Tosteson,  K.  D.  Paulsen,  "Automated  Region  Detection  Based 
Upon  Contrast  To  Noise  Ratio  in  Near-Infrared  Tomography,"  Applied  Optics  43(5)  1053-1062,  2004 

H.  Dehghani,  M.  M.  Doyley,  B.  W.  Pogue,  S.  Jiang,  J.  Geng  and  K.  D.  Paulsen,  "Breast  deformation  modeling  for  image 
reconstruction  in  near  infrared  optical  tomography"  Phys.  Med.  Biol.  Volume  49,  Number  7,  pp  1131-1146  2004. 

S.  Srinivasan,  B.  W.  Pogue,  H.  Dehghani,  S.  Jiang,  X.  Song  and  K.  D.  Paulsen,  "Improved  Quantification  of  Small  Objects  in  Near- 
Infrared  Diffuse  Optical  Tomography  JBO  9(6)  1161-1171  2004 

B.  W.  Pogue,  S.  Jiang,  H.  Dehghani,  C.  Kogel,  S.  Soho,  S.  Srinivasan,  X.  Song,  S.  P.  Poplack,  and  K.  D.  Paulsen,  "Characterization 
of  Hemoglobin,  Water  and  NIR  Scattering  in  Breast  Tissue:  Analysis  of  inter-subject  variability  and  menstrual  cycle  changes  relative 
to  lesions"  JBO,  9(3)  541-552,  2004. 

Ben  Brooksby,  Shudong  Jiang,  Christine  Kogel,  Marvin  Doyley,  Hamid  Dehghani,  John  B  Weaver,  Steven  P  Poplack,  Brian  W 
Pogue  and  Keith  D  Paulsen,  "Magnetic  Resonance  Guided  Near  Infrared  Tomography  of  the  Breast ",  Rev.  Sci.  Inst.  75  5262-5270 

2004 

Heng  Xu,  Brian  W.  Pogue,  Hamid  Dehghani,  Keith  D.  Paulsen,  Roger  Springett  and  Jeff  F.  Dunn,  "Absorption  and  scattering 
imaging  of  tissue  with  steady-state  second-differential  spectral- analysis  tomography"  Optics  letters  29(17)  2043-2045  2004 
Hamid  Dehghani,  Nirmal  Soni,  Ryan  Halter,  Alex  Hartov  and  Keith  D  Paulsen,  "Excitation  patterns  in  three-dimensional  electrical 
impedance  tomography"  Physiol.  Meas.  26  S185-S197  2005 

Heng  Xu,  Roger  Springett,  Hamid  Dehghani,  Brian  W  Pogue,  Keith  D  Paulsen  and  Jeff  F  Dunn,  "MRI  coupled  Broadband  Near 
Infrared  Tomography  System  for  Small  Animal  Brain  Studies ",  Applied  Optics  44(1 1)  2177-2188  2005 

Hamid  Dehghani,  Ben  Brooksby,  Brian  W  Pogue  and  Keith  D  Paulsen,  "Effect  of  Refractive  Index  on  Near  Infrared  Tomography  of 
the  Breast"  Applied  Optics  44(10)  1870-1878,  2005 

Subhadra  Srinivasan,  Brian  W  Pogue,  Shudong  Jiang,  Hamid  Dehghani  and  Keith  D  Paulsen,  "Spectrally  Constrained  Chromophore 
and  Scattering  NIR  Tomography  Improved  Quantification  and  Robustness  of  Reconstruction" ,  Applied  Optics  44(10)  1858-1869, 

2005 

Ben  Brooksby,  Shudong  Jiang,  Hamid  Dehghani,  Brian  W.  Pogue,  Keith  D.  Paulsen,  "Combining  near  infrared  tomography  and 
magnetic  resonance  imaging  to  study  in  vivo  breast  tissue:  implementation  of  a  Laplacian-type  regularization  to  incorporate  MR 
structure"  JBO  2005  Accepted 

Ben  Brooksby,  Subhadra  Srinivasan,  Shudong  Jiang,  Hamid  Dehghani,  Brian  W.  Pogue,  Keith  D.  Paulsen,  John  Weaver,  Christine 
Kogel  and  Steven  P.  Poplack,  "Spectral-prior  information  improves  Near-Infrared  diffuse  tomography  more  than  spatial-prior"  Optic 
Letters  2005  Accepted 

Subhadra  Srinivasan,  Brian  W.  Pogue,  Hamid  Dehghani,  Shudong  Jiang,  Xiaomei  Song,  Keith  D.  Paulsen,  "Effect  of  image 
reconstruction  bias  upon  spectroscopy-based  quantification  of  chromophores  in  nearinfrared  tomography"  In  Biomedical  Topical 
Meetings  on  CD-ROM  (The  Optical  Society  of  America,  Washington,  DC)  2004. 

Ben  A.  Brooksby,  Shudong  Jiang,  Gordon  Ehret,  Hamid  Dehghani,  Brian  W.  Pogue,  Keith  D.  Paulsen,  "Development  of  a  system  for 
simultaneous  MRI  and  Near-infrared  diffuse  tomography  to  diagnose  breast  cancer  In  Biomedical  Topical  Meetings  on  CD-ROM 
(The  Optical  Society  of  America,  Washington,  DC)  2004. 

Heng  Xu,  Hamid  Dehghani,  Brian  W.  Pogue,  Keith  D.  Paulsen,  Roger  Springett,  Jeff  F.  Dunn,  "Optical  tomography  system  based  on 
second-differential  spectroscopy  for  small  animal  brain  study  In  Biomedical  Topical  Meetings  on  CD-ROM  (The  Optical  Society  of 
America,  Washington,  DC)  2004. 

Brian  Pogue,  Shudong  Jiang,  Subhadra  Srinivasan,  Xiaomei  Song,  Hamid  Dehghani,  Keith  Paulsen,  Tor  Tosteson,  Christine  Kogel, 
Sandra  Soho,  Steven  P.  Poplack,  "Near-infrared  scattering  spectrum  differences  between  benign  and  malignant  breast  tumors 
measured  in  vivo  with  diffuse  tomography  In  Biomedical  Topical  Meetings  on  CD-ROM  (The  Optical  Society  of  America, 
Washington,  DC)  2004. 

Hamid  Dehghani,  Brian  W.  Pogue,  Marvin  M.  Doyley,  Jason  Geng,  Keith  D.  Paulsen,  "Breast  deformation  in  near  infrared  optical 
tomography  In  Biomedical  Topical  Meetings  on  CD-ROM  (The  Optical  Society  of  America,  Washington,  DC)  2004 
Ben  Brooksby,  Subhadra  Srinivasan,  Brian  W  Pogue,  Shudong  Jiang,  Hamid  Dehghani,  Christine  Kogel,  John  Weaver,  Steven 
Poplack,  Justin  D  Pearlman  and  Keith  D  Paulsen,  "Quantifying  adipose  and  fibroglandular  breast  tissue  properties  using  MRI- guided 
NIR  tomography ",  Proceedings  of  SPIE  5693,  2005 


8 


Phaneendra  K  Yalavarthy,  Hamid  Dehghani,  Brian  W  Pogue  and  Keith  D  Paulsen,  "Measurement  optimization  for  Near-Infrared 
optical  tomography ",  Proceedings  of  SPIE  5693,  2005 

Subhadra  Srinivasan,  Brian  W  Pogue,  Shudong  Jiang,  Hamid  Dehghani  and  Keith  D  Paulsen,  "Spectrally  Constrained  NIR 
tomography  for  breast  Imaging:  Simulations  and  Clinical  Results ",  Proceedings  of  SPIE  5693,  2005 

Appendices 

Four  papers  are  included  as  appendices,  which  are  as  a  direct  result  of  the  project  funding: 


9 


Institute  of  Physics  Publishing 


Physics  in  Medicine  and  Biology 


Phys.  Med.  Biol.  49  (2004)  1131-1145 


PII:  S003 1-9155(04)74363-8 


Breast  deformation  modelling  for  image 
reconstruction  in  near  infrared  optical  tomography 

Hamid  Dehghani1,  Marvin  M  Doyley2,  Brian  W  Pogue1,  Shudong  Jiang1, 
Jason  Geng3  and  Keith  D  Paulsen1 

1  Thayer  School  of  Engineering,  Dartmouth  College,  Hanover,  NH  03755,  USA 

2  Department  of  Radiology,  Dartmouth  Medical  School,  Hanover,  NH  03755,  USA 

3  Genex  Technologies  Inc.,  10605  Concord  Street,  Suite  500  Kensington,  MD  20895-2504,  USA 

Received  8  October  2003 
Published  18  March  2004 

Online  at  stacks.iop.org/PMB/49/1131  (DOI:  10.1088/0031-9155/49/7/004) 

Abstract 

Near  infrared  tomography  (NIR)  is  a  novel  imaging  technique  that  can  be 
used  to  reconstruct  tissue  optical  properties  from  measurements  of  light 
propagation  through  tissue.  More  specifically  NIR  measurements  over  a  range 
of  wavelengths  can  be  used  to  obtain  internal  images  of  physiologic  parameters 
and  these  images  can  be  used  to  detect  and  characterize  breast  tumour.  To 
obtain  good  NIR  measurements,  it  is  essential  to  have  good  contact  between 
the  optical  fibres  and  the  breast  which  in-turn  results  in  the  deformation  of  the 
breast  due  to  the  soft  plasticity  of  the  tissue.  In  this  work,  a  tissue  deformation 
model  of  the  female  breast  is  presented  that  will  account  for  the  altered  shape 
of  the  breast  during  clinical  NIR  measurements.  Using  a  deformed  model  of 
a  breast,  simulated  NIR  data  were  generated  and  used  to  reconstruct  images 
of  tissue  absorption  and  reduced  scatter  using  several  assumptions  about  the 
imaging  domain.  Using  either  a  circular  or  irregular  2D  geometry  for  image 
reconstruction  produces  good  localization  of  the  absorbing  anomaly,  but  it 
leads  to  degradation  of  the  image  quality.  By  modifying  the  assumptions  about 
the  imaging  domain  to  a  3D  conical  model,  with  the  correct  diameter  at  the 
plane  of  NIR  measurement,  significantly  improves  the  quality  of  reconstructed 
images  and  helps  reduce  image  artefacts.  Finally,  assuming  a  non-deformed 
breast  shape  for  image  reconstruction  is  shown  to  lead  to  poor  quality  images 
since  the  geometry  of  the  breast  is  greatly  altered,  whereas  using  the  correct 
deformed  geometry  produces  the  best  images. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Near  infrared  (NIR)  tomography  is  an  emerging  alternative  imaging  method  used  to 
image  physiologic  parameters  of  biological  tissue  in  vivo  such  as  water  and  haemoglobin. 
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Measurements  of  light  propagation  (600-900  nm)  within  tissue  can  be  used  to  map  internal 
chromophore  concentrations  within  tissue.  Light  is  transmitted  through  tissue  using  multiple 
input  and  output  locations  on  the  surface  of  the  region  to  be  imaged,  similar  to  a  fan- 
beam  x-ray  computed  tomography  geometry,  but  using  optical  fibres  for  delivery  and  pick 
up  of  the  light  signals.  The  intensity  and  path-length  distributions  of  the  exiting  photons 
provide  information  about  the  optical  properties  of  the  transilluminated  tissue  using  a  model- 
based  interpretation  where  photon  propagation  is  simulated  by  diffusion  theory.  Using  these 
techniques  implemented  in  hardware  and  software,  NIR  optical  tomography  becomes  an 
inherently  three-dimensional  imaging  method  and  is  used  to  reconstruct  physiologically 
relevant  chromophore  distributions  from  the  region  under  investigation  (Eda  et  al  1999, 
Fantini  et  al  1999,  Boas  et  al  2001,  Hebden  et  al  2001,  Pogue  et  al  2001,  McBride  et  al  2002, 
Dehghani  et  al  2003a,  2003b,  2003c).  The  main  interest  in  this  study  lies  in  the  ability  to 
detect  and  characterize  tumours  within  the  female  breast  (Pogue  et  al  2001,  Brooksby  et  al 
2003).  Since  the  absorption  and  scattering  of  light  in  tissue  is  a  function  of  its  optical 
properties,  and  hence  its  physiological  state,  the  aim  is  to  obtain  images  of  internal  optical 
absorption  coefficient  /xa,  reduced  scattering  coefficient  /xs’  and  ultimately  images  of  total 
haemoglobin  and  oxygen  saturation  distributions.  These  images  should  in  principle  provide 
information  about  the  physiological  state  of  the  tissue  under  investigation  and  help  identify 
and  characterize  tumours  within  the  breast.  However,  in  recent  years  it  has  become  apparent 
that  the  external  shape  of  the  tissue  can  dramatically  impact  the  quality  of  the  reconstruction, 
and  this  is  largely  due  to  mismatch  between  the  model  prediction  and  the  actual  shape  of  the 
tissue  boundary.  In  this  study,  a  focused  effort  has  been  put  forward  to  examine  the  extent  to 
which  external  boundary  changes  will  degrade  image  quality  and,  more  importantly,  how  this 
could  be  corrected  with  appropriate  model-based  analysis. 

To  obtain  clinical  measurements  with  a  sufficient  signal-to-noise  ratio,  it  is  key  to  ensure 
that  good  contact  exists  between  the  optical  fibres  and  the  tissue.  More  specifically  to  our 
studies,  the  patient  lies  prone  on  the  measurement  bed,  which  contains  a  single  opening  for 
the  breast.  The  breast  is  suspended  freely  through  this  opening,  below  which  the  optical  fibres 
are  brought  into  contact  with  the  breast.  Normally,  the  optical  fibre  arrangement  consists  of 
a  total  of  48  fibres,  16  fibres  in  three  separate  planes.  However,  for  the  purpose  of  this  work 
only  one  plane  of  16  fibres  is  considered.  The  fibres  need  to  make  full  contact  with  the  breast 
for  a  good  and  adequate  NIR  measurement. 

The  breast  is  a  soft  tissue,  which  will  deform  and  alter  its  shape  on  the  application  of 
external  pressure.  The  amount  of  deformation  is  a  function  of  the  tissue  mechanical  properties 
and  the  amount  of  displacement  and  the  external  pressure  applied  by  the  optical  fibre  array.  In 
most  image  reconstruction  algorithms,  the  general  assumption  is  made  that  the  region  under 
investigation  is  a  uniform  circular  (2D)  or  conical  or  cylindrical  (3D)  domain.  Little  work  has 
been  done  to  evaluate  the  effect  of  any  incorrect  geometry  in  image  reconstruction.  For  simple 
symmetric  geometries,  for  example  a  3D  cylinder,  work  has  been  described  to  accurately 
calibrate  the  data  to  account  for  2D/3D  mismatch  (Hillman  et  al  2000).  More  recently,  work 
has  been  done  with  regard  to  neonatal  head  imaging  that  a  slight  change  to  the  geometry 
model,  both  in  terms  of  geometry  and  mesh,  will  result  into  large  change  in  modelled  data 
(Gibson  et  al  2003a). 

In  this  work,  the  effect  of  such  assumptions  is  investigated  by  creating  a  deformation  of 
the  breast  model  in  3D  to  generate  simulated  clinical  data.  Images  are  then  reconstructed  using 
various  assumptions  regarding  the  imaging  domain  including  a  circular  or  irregular  2D  models 
as  well  as  a  3D  conical  shaped  model  and  the  non-deformed  breast  model  to  evaluate  image 
quality.  An  important  part  of  this  modelling  effort  is  the  recent  advances  which  have  been  made 
in  modelling  soft  tissue  deformation  (Paulsen  et  al  1999,  Doyley  et  al  2000,  Van  Houten  et  al 
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2000).  Beyond  this,  it  is  also  becoming  established  that  soft  tissue  elastic  properties  can  be 
directly  measured  with  elastography  imaging  in  vivo ,  thereby  providing  the  key  information 
needed  for  predicting  the  deformation  of  tissue  under  force  or  displacement  conditions.  In  this 
study,  this  displacement  modelling  was  used  along  with  the  NIR  tomography  modelling 
to  create  a  comprehensive  three-dimensional  prediction  of  how  to  approach  appropriate 
modelling  of  deformed  tissue  that  is  being  imaged  with  NIR  tomography. 

2.  Theory 

2.1.  Breast  deformation  model 

Soft  tissues  exhibit  nonlinear  elastic  behaviour  (Fung  1993).  Nevertheless,  it  can  be  considered 
a  linear  elastic  material  in  situations  where  the  deforming  forces  produce  infinitesimal 
deformations  (i.e.  ^5%)  (Timoshenko  and  Goodier  1970).  For  the  purpose  of  this  study, 
the  breast  was  modelled  as  linear  isotropic  pseudo-incompressible  (i.e.  Poisson’s  ratio  (v)  = 
0.495  (Fung  1993)).  Under  these  assumptions  and  ignoring  internal  body  force,  the  governing 
elasticity  equations  for  quasi-static  deformation  are  given  by  Timoshenko  and  Goodier  (1970) 
and  Fung  (1993) 


(A  +  /x)  V  (V  •  u )  +  /xV  u  —  0 
for  internal  nodes  in  domain  Q,  and 

((A  +  /x)V(V  •  u)  +  /xV2w)  •  h  =  h 


(1) 


(2) 


for  nodes  on  the  boundary  8Q. 

Here  h  represents  a  unit  vector  directed  outwards  from  Q,  and  h  represents  the  traction  on 
the  surface  or  boundary  of  the  breast.  Note  that  u  represents  the  displacement  components  in 
all  coordinate  directions,  and  /x  and  y  are  Lame’s  elastic  constants.  For  an  isotropic  medium 
these  constants  are  related  to  the  more  familiar  Young’s  modulus  (E)  and  the  Poisson’s  ratio 
(v)by 

E  vE 

/x  =  -  A  =  - .  (3) 

^  2(1 +  v)  (1  +  v)(l  —  2v) 

The  first  Lame’s  elastic  constant  (i.e.  /x)  is  generally  known  as  shear  modulus. 

It  should  be  stated  that  in  this  study,  the  breast  was  assumed  to  be  traction  free  (i.e.  no 
other  forces  are  associated)  and  internal  body  forces  were  neglected,  thus  the  problem  is  solved 
by  imposing  a  prescribed  displacement  (i.e.  induced  when  the  optical  fibres  are  coupled  to  the 
breast)  as  described  in  Doyley  et  al  (1999). 


2.2.  Light  propagation  model 

Under  the  assumption  that  scattering  dominates  absorption  for  NIR  light  in  tissue,  the 
Boltzmann  transport  equation  can  be  simplified  to  the  diffusion  approximation,  which  in 
the  frequency  domain  is  given  by 


/  i<x>\ 

—V  •  DVO(r,  co)  +  ypLa  +  — J  d>(r,  co)  =  qo(r,  co) 


(4) 


where  qo(r ,  co)  is  an  isotropic  source,  0(r,  co)  is  the  photon  fluence  rate  at  position  r  and 
D  =  1  is  the  diffusion  coefficient.  We  use  the  Robin  (Type  III)  boundary  condition: 

OQ/)  +  -fi-  VO(y)  =  0 
a 


(5) 
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where  a  is  a  term  which  incorporates  reflection  as  a  result  of  refractive  index  mismatch 
(Schweiger  et  al  1995,  Dehghani  et  al  2003c)  at  the  boundary,  and  h  is  the  outward  pointing 
normal  to  the  boundary  (8Q)  at  y. 

We  assume  that  the  data  are  represented  by  a  nonlinear  operator  y*  =  F[/jia,  D],  where 
our  data  y*  are  a  complex  vector  having  a  real  and  imaginary  components,  which  are  mapped 
to  log  amplitude  and  phase  shift  in  measurement.  Then  the  image  reconstruction  method  is 
posed  as  a  solution  to  the  following  expression: 

(Afl,  D)  =  argmin||(y*  -  F(fia,  D))\\  (6) 

M a,K 

where  ||-||  is  the  weighted  L2-norm,  representing  the  square  root  of  the  sum  of  the  squared 
elements,  jla  is  a  vector  of  the  absorption  coefficients  and  D  is  a  vector  of  the  diffusion 
coefficients.  The  magnitude  of  this  is  sometimes  referred  to  as  the  projection  error  and 
provides  a  value  for  determining  the  convergence  of  the  iterative  reconstruction  algorithm. 

In  this  study,  a  finite-element  method  (FEM)  is  used  as  a  general  and  flexible  method  for 
solving  the  forward  problem  in  arbitrary  geometries  (Arridge  et  al  1993,  Jiang  et  al  1996).  In 
the  inverse  problem,  where  the  goal  is  to  recover  internal  optical  property  distributions  from 
boundary  measurements,  it  is  assumed  that  tia(r)  and  D( r)  are  expressed  in  a  basis  with  a 
limited  number  of  dimensions  (less  than  the  dimension  of  the  finite  element  system  matrices). 
A  number  of  different  strategies  for  defining  reconstruction  bases  are  possible;  in  this  paper  a 
linear  pixel  basis  (Schweiger  and  Arridge  1999)  is  used.  To  find  (pLa,  D)  in  equation  (6)  we 
have  used  a  Levenberg-Marquardt  algorithm,  where  we  repeatedly  solve: 

a  =  JT(JJT  +  pl)~lb  (7) 

where  b  is  the  data  vector,  b  =  [  y*  —  F\ii„,  D]]r,  a  is  the  solution  update  vector, 
a  =  [8D;  8fla].  Here,  p  is  the  regularization  factor  and  /  is  the  Jacobian  matrix  for  our 
model  which  is  calculated  using  the  Adjoint  method  (Arridge  and  Schweiger  1995,  Dehghani 
et  al  2003b,  2003c). 

3.  Methods 

3.1.  Breast  deformation 

A  volume  mesh  of  a  female  breast  of  a  volunteer  was  created  from  surface  image  data  that  was 
acquired  using  a  3D  surface  camera  (Rainbow  3D  Camera,  Genex  Technologies,  Kensington 
MD).  The  3D  camera  projects  structural  illumination  patterns  onto  the  object  and  calculates 
3D  surface  profile  described  by  over  300  000  data  points  (Geng  1996,  Galdino  et  al  2002).  A 
volume  mesh  is  then  generated  using  the  Delaunay  algorithm.  The  mesh  has  a  geometry  of 
130  x  136  x  60  mm,  figure  1(a),  and  contained  15  501  nodes  corresponding  to  61  171  linear 
tetrahedral  elements.  The  diameter  of  the  breast  at  its  mid-plane  where  optical  fibre  array  will 
be  applied  is  approximately  88  mm.  To  calculate  the  deformation  due  to  16  equally  spaced 
optical  fibres  being  applied  at  the  mid-plane  of  the  breast,  i.e.  z  =  —  30  mm,  it  is  assumed  that 
each  optical  fibre  pushed  the  breast  so  that  the  final  breast  diameter  at  z  =  —  30  mm  is  70  mm, 
and  that  the  diameter  of  each  optical  fibre  is  6  mm.  The  modelled  elastic  properties  of  tissue, 
equation  (3),  were  assumed  as  isotropic  and  homogenous  with  Young’s  modulus  of  20  kPa 
(Krouskop  et  al  1998)  and  Poisson’s  ratio  of  0.495  (Fung  1993).  Further,  it  was  assumed  that 
the  topmost  part  of  the  mesh,  i.e.  z  =  0  mm  was  not  allowed  to  move  since  it  is  connected 
to  the  chest.  Using  this  applied  displacement  as  a  boundary  condition,  the  displacement  at 
all  nodes  due  to  the  application  of  the  optical  fibres  was  calculated  and  a  deformed  mesh  was 
created,  figure  1(b). 
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Figure  1.  Volume  mesh  of  the  (a)  normal  suspended  breast  and  (b)  the  deformed  mesh  after  the 
application  of  the  optical  fibre  array. 


3.2.  Simulation  of  data  from  deformed  breast 

In  order  to  accurately  simulate  the  clinical  settings,  two  localized  anomaly  regions  were  placed 
within  the  mid-plane  of  the  normal  breast,  both  at  z  =  —30  mm,  figure  2(a).  First  was  an 
absorbing  anomaly,  19  mm  from  the  surface  and  a  radius  10  mm  with  a  coefficient  value  of 
0.03  mm-1.  Second  was  a  reduced  scattering  anomaly,  also  19  mm  from  the  surface  but  with 
a  radius  of  7.5  mm  and  a  value  of  3  mm-1.  These  anomalies  were  chosen  as  they  represent 
the  size  and  contrast  we  aim  to  image  and  so  that  absorption  and  scatter  separation  can  also  be 
investigated.  The  background  optical  properties  were  modelled  with  an  absorption  coefficient 
of  0.01  mm-1  and  a  reduced  scatter  coefficient  of  1.0  mm-1.  For  the  deformed  breast  model 
the  same  optical  properties  of  the  normal  breast  were  assumed  (the  anomalies  were  also 
assumed  to  have  the  same  elastic  properties  of  the  background).  Once  the  deformation  of  the 
model  was  calculated,  it  was  assumed  that  the  anomalies  were  free  to  move,  depending  on 
the  elastic  properties  of  the  breast,  figure  2(b).  It  is  interesting  to  note  that  the  displacement  of 
the  anomalies  is  not  so  much  dependent  on  the  mesh  density  of  the  breast,  but  more  dependent 
on  the  mechanical  properties,  applied  pressure  and  non-symmetric  nature  of  the  breast. 

Using  the  deformed  mesh,  together  with  the  displaced  anomalies,  data  were  simulated  for 
16  equally  circularly  spaced  optical  fibres  placed  at  z  =  —30  mm.  Amplitude  and  phase  data 
were  simulated  at  100  MHz,  and  1%  noise  was  added.  These  data  were  then  used  as  simulated 
patient  data  in  the  following  sections.  In  addition,  to  allow  data  calibration  as  done  using 
clinical  data,  and  discussed  elsewhere  (Dehghani  et  al  2003c,  McBride  et  al  2003),  the  data 
were  simulated  for  16  equally  spaced  optical  fibres  placed  in  a  circle  around  the  mid-plane  of  a 
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z  =  -50  mm  z  =  -40  mm  z  =  -30  mm  z  =  -20  mm  z  =  -10  mm 


Absorption 


z  =  -50  mm  z  =  -40  mm  z  =  -30  mm  z  =  -20  mm  z  =  -10  mm 


z  =  -50  mm  z  =  -40  mm  z  =  -30  mm  z  =  -20  mm  z  =  -10  mm 


Absorption 


z  =  -50  mm  z  =  -40  mm  z  =  -30  mm  z  =  -20  mm  z  =  -10  mm 


Figure  2.  2D  coronal  slices  through  (a)  the  normal  breast  mesh  and  (b)  the  deformed  breast  mesh, 
showing  the  position  of  the  anomalies.  The  most  right-hand  slice  is  near  the  chest  while  the  most 
left-hand  slice  is  near  the  nipple. 
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Figure  3.  2D  simultaneous  reconstruction  of  absorption  and  reduced  scatter  from  deformed  breast 
model  simulated  data.  The  mesh  used  for  the  reconstruction  is  a  circular  mesh  whose  diameter  is 
the  same  as  the  optical  fibre  array  diameter  used  to  deform  and  simulate  data  NIR  from  the  breast 
model. 


non-deformable  cylindrical  model  with  a  known  homogenous  background  optical  properties. 
The  same  calibration  data  file  has  been  used  for  all  presented  reconstructions. 

3.3.  Image  reconstruction  using  2D  meshes 

Using  the  computed  data  for  the  deformed  mesh  and  the  reference  phantom,  as  simulated 
measurements  with  1%  random  noise  added,  the  data  were  calibrated  (Dehghani  et  al  2003c, 
McBride  et  al  2003)  and  images  were  reconstructed  using  two  different  2D  meshes.  Assuming 
that  the  imaging  domain  was  a  circular  domain  with  a  radius  equal  to  the  radius  of  the  optical 
fibre  array,  i.e.  35  mm,  images  were  reconstructed  and  the  results  are  shown  in  figure  3.  The 
circular  mesh  used  contained  1785  nodes  corresponding  to  3418  linear  triangular  elements. 
For  the  reconstruction  basis  (basis  used  for  the  update  of  the  optical  parameters),  a  20  x  20 
regular  grid  (piecewise  linear)  was  used  and  the  initial  regularization  parameter  was  set  to  100. 
Images  shown  here,  and  all  following  sections  are  at  iteration  level  where  the  projection  error 
was  within  5%  of  the  previous  iteration.  In  the  2D  reconstruction  cases,  this  projection  error 
was  reached  typically  in  the  7th  iteration. 

It  is  evident  from  figure  1(b)  that  the  breast  at  the  plane  of  measurement,  i.e.  z  =  —  30  mm, 
was  no  longer  circular.  In  order  to  evaluate  if  a  more  correct  a  priori  information  regarding 
the  2D  boundary  at  the  plane  of  measurement  would  be  useful,  a  2D  irregular  mesh,  with  a 
boundary  whose  profile  is  the  same  as  the  boundary  at  z  =  —30  mm  for  the  deformed  mesh, 
was  created  and  used  for  image  reconstruction.  The  resulting  images  reconstructed  on  this 
mesh  are  shown  in  figure  4.  The  irregular  mesh  used  contained  2405  nodes  corresponding 
to  4619  linear  triangular  elements.  For  the  reconstruction  basis,  a  20  x  20  regular  grid  was 
used  and  the  initial  regularization  parameter  was  set  to  100.  Thus  in  this  case,  the  irregular 
boundary  was  matched,  but  still  the  image  reconstruction  was  approximated  by  a  2D  forward 
solution. 

3.4.  Image  reconstruction  using  3D  conical  mesh 

In  the  next  section,  we  examined  the  improvement  if  a  true  3D  forward  calculation  was 
assumed,  but  with  a  regular  geometry.  In  a  previous  work,  it  was  hypothesized  that  given 
a  set  of  3D  patient  data,  those  images  could  be  reconstructed  with  reasonable  accuracy  if 
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Figure  4.  2D  simultaneous  reconstruction  of  absorption  and  reduced  scatter  from  deformed  breast 
model  simulated  data.  The  mesh  used  for  the  reconstruction  is  an  irregular  whose  boundary  is  the 
same  as  the  boundary  of  the  deformed  breast  mesh  at  z  =  30  mm. 


Figure  5.  The  conical  shaped  volume  mesh.  The  geometry  of  the  mesh  is  such  that  the  diameter  at 
the  mid-plane  corresponds  to  the  diameter  of  the  optical  fibre  array  used  to  deform  and  simulate  NIR 
data  of  the  breast  model.  The  diameter  at  20  mm  above  and  20  mm  below  the  NIR  measurement 
plane  of  the  deformed  breast  were  also  used  to  set  the  upper  and  lower  bound  of  the  cone  model. 


the  breast  was  assumed  to  be  a  conical  shaped  model  (Dehghani  et  al  2003c).  To  create  a 
conical  mesh,  the  information  about  the  diameter  of  the  measurement  fibre  array  as  well  as 
the  diameter  of  the  breast  at  20  mm  above  and  20  mm  below  the  measurement  plane  was 
used,  which  would  include  all  boundary  information  other  than  effects  due  to  tissue  bulging. 
Following  this  procedure,  a  conical  shaped  volume  mesh  was  created,  as  shown  in  figure  5, 
which  shows  the  upper  chest  wall  and  the  confined  partial-conical  breast  shape  as  a  regular 
geometry  approximation  to  the  circularly  compressed  breast  shape.  The  mesh  contained  9332 
nodes  corresponding  to  42  28 1  linear  tetrahedral  elements  and  was  created  using  NETGEN 
(Schoberl).  Using  this  mesh  and  the  simulated  data  for  the  deformed  breast  mesh,  images 
were  reconstructed  and  the  results  are  shown  in  figure  6. 

For  the  reconstruction  basis,  a  20  x  20  x  10  regular  grid  was  used  and  the  initial 
regularization  parameter  was  set  to  100.  The  images  shown  here  and  all  through  the  following 
section  are  at  iteration  level  where  the  projection  error  was  within  5%  of  the  previous  iteration. 
These  presented  from  the  3D  cases  are  from  the  14th  iteration. 

3.5.  Image  reconstruction  using  the  normal  and  deformed  breast  mesh 

Assuming  that  correct  information  regarding  the  breast  is  available,  before  the  application 
of  the  optical  fibre  array,  the  normal,  non-deformed  mesh  was  used  as  shown  in  figure  1(a), 
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Reduced  Scatter 

Figure  6.  3D  simultaneous  reconstruction  of  absorption  and  reduced  scatter  from  deformed  breast 
model  simulated  data.  The  mesh  used  for  the  reconstruction  is  the  mesh  shown  in  figure  5.  The 
most  right-hand  slice  is  near  the  chest  while  the  most  left-hand  slice  is  near  the  nipple. 


as  the  reconstruction  mesh  for  the  images  shown  in  figure  7.  The  reconstruction  basis  and 
parameters  used  were  the  same  as  for  the  conical  model  already  described.  Finally,  to  show 
the  best  possible  reconstruction,  the  deformed  mesh,  figure  1(b)  was  used  for  reconstruction 
of  images  and  these  are  shown  in  figure  8. 


4.  Results 

The  deformed  breast  mesh,  as  a  result  of  applying  a  circular  optical  fibre  array  at  z  =  30  mm, 
is  shown  in  figure  1(b).  It  is  evident  that  the  mesh  has  been  compressed  at  the  plane  of  applied 
displacement  with  slight  bulging  in-between  each  optical  fibre.  The  breast  mesh  has  also 
extended  in  the  z  direction  near  the  nipple  to  maintain  a  constant  volume  as  the  nodes  on  the 
chest  wall  are  assumed  fixed.  Although  the  mechanical  properties  of  the  tissue  are  assumed 
constant,  it  is  interesting  to  note  that  the  modelled  anomalies  have  also  been  displaced  due  to 
the  deformation,  figures  2(a)  and  (b).  Both  the  anomalies  have  been  reduced  in  diameter,  but 
have  extended  in  the  z  direction  as  a  result  of  the  applied  deformation. 

Two-dimensional  images  reconstructed  using  the  simulated  deformed  breast  data  are 
shown  in  figures  3  and  4.  In  both  cases  where  a  circular  or  an  irregular  boundary  was 
used,  good  images  were  recovered  in  terms  of  localization  of  the  absorption  anomaly  (within 
1.2  mm  for  the  circular  model,  and  2.0  mm  for  the  irregular  model).  The  recovery  of  the 
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Figure  7.  3D  simultaneous  reconstruction  of  absorption  and  reduced  scatter  from  deformed  breast 
model  simulated  data.  The  mesh  used  for  the  reconstruction  is  the  mesh  shown  in  figure  1(a).  The 
most  right-hand  slice  is  near  the  chest  while  the  most  left-hand  slice  is  near  the  nipple. 


reduced  scattering  anomaly  was  not  as  good  as  the  absorption,  and  both  sets  of  images  contain 
large  boundary  artefacts.  The  background  values  of  both  the  absorption  and  reduced  scatter 
were  close  to  the  expected  value,  however,  the  target  absorption  value  is  less  than  50%  than 
expected. 

Three-dimensional  images  reconstructed  assuming  a  conical  shaped  breast  are  shown  in 
figure  6.  In  that  case,  both  the  absorption  and  the  scattering  anomalies  were  recovered  with 
good  localization  (within  2.00  mm  for  both  the  absorption  and  scattering  objects)  and  with 
little  boundary  artefacts.  Since  the  measured  data  were  about  the  mid-plane  of  the  breast, 
and  the  modelled  cone  geometry  also  had  the  optical  fibres  around  its  mid-plane,  the  objects 
recovered  are  centred  at  z  =  —30  mm.  Here  the  background  values  of  absorption  and  scatter 
were  a  little  lower  than  expected  values,  however,  the  absorption  anomaly  had  a  maximum 
value  of  0.014  mm-1  and  the  scattering  object  had  a  maximum  value  of  2  mm-1. 

Images  reconstructed  using  the  assumption  of  correct  non-deformed  breast  geometry  are 
shown  in  figure  7.  In  this  case,  although  the  absorbing  anomaly  has  been  reconstructed,  the 
reduced  scattering  images  contain  large  artefacts.  The  background  absorption  and  scatter 
values  are  about  0.006  mm-1  and  0.5  mm-1,  respectively,  while  the  maximum  absorption  for 
the  anomaly  is  0.014  mm-1.  Finally,  images  reconstructed  assuming  correct  3D  deformed 
boundary  are  shown  in  figure  8.  Here,  both  the  absorbing  and  scattering  regions  have  been 
recovered  within  less  than  1.00  mm  of  expected  location,  giving  great  localization  accuracy. 
There  exists  some  cross  talk  between  the  absorbing  and  scattering  anomalies.  The  background 
absorption  and  scatter  values  are  about  0.01  mm-1  and  0.98  mm-1,  respectively,  while  the 
maximum  absorption  and  scatter  values  for  the  anomaly  is  0.013  mm-1  and  1.3  mm-1, 
respectively. 
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Absorption 


z  =  -50  mm  z  =  -40  mm  z  =  -30  mm  z  =  -20  mm  z  =  -10  mm 


Reduced  Scatter 


Figure  8.  3D  simultaneous  reconstruction  of  absorption  and  reduced  scatter  from  deformed  breast 
model  simulated  data.  The  mesh  used  for  the  reconstruction  is  the  mesh  shown  in  figure  1(b).  The 
most  right-hand  slice  is  near  the  chest  while  the  most  left-hand  slice  is  near  the  nipple. 


5.  Discussions 

The  deformation  of  the  breast  due  to  the  application  of  16  circular  equally  spaced  optical 
fibres  has  been  modelled  using  a  computational  mechanic  model.  Using  this  deformed 
model  which  also  contains  a  single  absorbing  and  a  single  scattering  perturbations,  data  were 
simulated  and  1  %  noise  was  added  to  more  adequately  represent  clinical  measurements.  Using 
these  simulated  data,  reconstructed  images  of  absorption  and  reduced  scatter  were  generated 
simultaneously  using  various  different  geometries  of  the  mesh  upon  which  the  reconstruction 
would  take  place.  In  the  first  case,  it  was  assumed  that  the  model  used  for  reconstruction  was 
circular  with  a  radius  equal  to  the  radius  of  the  optical  fibre  array  after  breast  compression, 
which  was  35  mm.  The  reconstructed  images  have  shown  good  accuracy  in  the  localization  of 
the  absorbing  object  with  a  maximum  absorption  value  of  0.014  mm-1.  The  recovery  of  the 
scattering  objects  is  much  poorer  than  the  absorbing  anomaly,  which  can  be  attributed,  perhaps, 
to  its  small  size  (15  mm  diameter).  There  exist  boundary  artefacts  in  both  the  absorbing  and 
scattering  images.  It  is  possible  to  improve  the  quality  of  the  reconstructed  2D  images  by 
using  several  methods,  including  the  use  of  various  regularization  parameters  and/or  use  of 
a  priori  information.  A  good  method  to  improve  the  reconstructed  images  in  this  case  would 
be  to  segment  the  reconstructed  images  shown  in  figure  3,  and  use  these  as  a  priori  information 
to  more  accurately  obtain  quantitative  results.  The  application  of  such  algorithm  is  discussed 
elsewhere,  and  has  been  found  to  provide  superior  results  (Srinivasan  et  al  2004).  In  the  case 
of  knowing  and  using  correct  2D  boundary  information  in  image  reconstruction,  figure  4,  no 
useful  improvement  was  seen.  Although  correct  information  was  added  to  the  reconstruction 
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with  regard  to  the  boundary,  the  image  reconstruction  algorithm  was  still  constrained  to 
be  2D.  While  it  has  been  shown  many  times  that  the  2D  geometry  can  be  used  to  recover 
accurate  image  data  from  NIR  tomography  with  a  circular  boundary  (Pogue  and  McBride 
1999),  it  is  not  surprising,  perhaps,  to  find  that  this  simplification  does  not  continue  to  work 
as  the  exterior  boundary  becomes  more  complex  in  shape.  The  reconstruction  mesh  is  no 
longer  symmetric  and  contains  many  irregular  edges,  which  gives  rise  to  high  frequency 
noise  at  the  boundary,  and  those  changes  occur  symmetrically  around  each  optical  fibre 
making  the  calculation  even  more  dependent  upon  the  3D  geometry  out  of  the  plane  of 
imaging. 

Next  for  the  case  where  the  imaging  domain  was  assumed  to  be  conical  in  shape,  the 
reconstructed  images  were  shown  in  figure  6.  Here  it  was  assumed  that  the  mid-plane  of 
the  cone  had  a  radius  equal  to  the  radius  of  the  optical  fibre  array  after  breast  compression. 
Further,  it  was  assumed  that  the  measurements  of  the  diameter  of  the  breast  20  mm  above 
and  20  mm  below  the  measurement  plane  were  known.  Using  these  data,  a  conical  shaped 
mesh  could  be  created  for  any  patient  exam  as  shown  in  figure  5.  The  reconstructed  images 
had  recovered  both  the  absorbing  and  scattering  anomalies  with  reasonable  success.  The 
calculated  background  values  for  absorption  and  scatter  were  0.007  mm-1  and  0.94  mm-1, 
respectively.  There  exists  a  small  cross  talk  between  the  absorption  and  reduced  scatter  images, 
plus  artefacts  at  the  boundary  of  the  images.  The  qualitative  and  quantitative  accuracy  of  the 
conical  shaped  model  reconstruction  are  better  than  the  2D  circular  case,  simply  because  a 
more  accurate  3D  model,  rather  than  a  2D  model  is  used. 

In  the  case  where  it  is  assumed  that  the  exact  breast  geometry  was  known,  but  have 
ignored  information  regarding  the  compression  due  to  the  optical  fibres,  the  resulting  images 
are  shown  in  figure  7.  Here,  although  the  absorption  anomaly  was  recovered  with  relatively 
good  accuracy,  the  reduced  scattering  image  contains  a  large  artefact.  Furthermore,  the 
calculated  background  values  for  absorption  and  scatter  are  0.006  mm-1  and  0.5  mm-1, 
respectively,  which  are  much  lower  than  expected.  These  images  are,  perhaps,  not  as  accurate 
and  useful  as  the  conical  geometry  case,  since,  although  there  is  a  more  accurate  3D  model, 
the  diameter  of  the  breast  was  not  correct,  particularly  within  the  measurement  plane,  i.e.  z  = 
—30  mm.  Finally,  for  the  case  of  using  a  geometrically  correct  deformed  mesh  to  reconstruct, 
the  images  are  shown  in  figure  8.  Here,  the  reconstructed  images  of  absorption  and  scatter 
were  found  with  superior  localization  accuracy.  The  scattering  object  has  also  been  recovered 
in  this  case.  However,  the  quantitative  accuracy  of  the  recovered  anomalies  is  not  as  good 
with  maximum  values  for  absorption  and  scattering  objects  of  0.013  mm-1  and  1.3  mm-1, 
respectively.  This  low  quantitative  accuracy  has  been  reported  elsewhere  and  is  a  common 
problem  in  3D  NIR  imaging  algorithms  (Dehghani  et  al  2003b,  Gibson  et  al  2003b),  which 
may  be  solved  using  of  multistage  reconstruction  algorithms  (Srinivasan  et  al  2004)  or  with 
inclusion  of  a  priori  data  (Brooksby  et  al  2003).  Finally,  using  the  correct  model  with  correct 
information  regarding  the  deformed  boundary  has  produced  images  with  very  little  or  no 
artefacts. 


6.  Conclusions 

In  this  work,  we  have  presented  a  breast  deformation  model  that  would  account  for  the 
change  in  shape  and  geometry  of  the  breast  due  to  the  application  of  a  circular  array  of 
optical  fibres  for  NIR  measurements.  The  proposed  model,  at  present,  assumes  that  the  breast 
has  homogenous  mechanical  properties,  which,  although  not  accurate,  provide  a  good  initial 
estimate  for  modelling  of  any  deformation.  Future  studies  may  focus  on  using  spatially 
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distributed  mechanical  properties,  as  these  can  be  imaged  with  good  spatial  resolution  with 
MR  elastography  (Paulsen  et  al  1999,  Van  Houten  et  al  2000,  Weaver  et  al  2001,  Doyley  et  al 
2003).  Also,  in  this  work,  we  have  deformed  the  breast  more  than  adequately  needed  by 
19  mm,  whereas  in  a  real  clinical  setting,  the  deformation  will  be  of  lower  magnitude,  typically 
about  5-10  mm.  In  practice  for  our  imaging  exams,  we  currently  attempt  to  deform  the  breast 
as  minimally  as  possible,  however,  there  are  strategic  benefits  to  compressing  in  that  the  signal 
transmitted  can  be  higher,  and  there  can  be  pressure-induced  changes  which  might  provide 
meaningful  contrast  about  what  the  tissue  is  composed  of  (Jiang  et  al  2003).  Nevertheless  we 
have  chosen  such  a  large  deformation  to  provide  us  with  a  worst-case  scenario  for  NIR  image 
reconstruction,  and  to  examine  if  it  is  feasible  to  exploit  larger  magnitude  deformations 
in  clinical  studies  while  not  compromising  the  integrity  of  the  data  collected  from  the 
breast. 

Using  simulated  data  from  the  deformed  breast  that  also  contains  optical  anomalies,  we 
have  reconstructed  images  assuming  various  reconstruction  geometries.  In  the  first  case 
of  assuming  that  the  reconstruction  mesh  is  either  a  circular  or  irregular  2D  mesh,  the 
reconstructed  images  have  shown  good  accuracy  in  recovering  the  location  of  the  absorbing 
anomaly.  The  2D  reconstructed  scattering  images  are  not  as  good  as  expected,  which  is  perhaps 
due  to  the  fact  that  the  scattering  object  had  a  small  diameter  of  15  mm  and  a  contrast  which 
was  three  times  the  background,  and  thus  the  image  being  dominated  by  boundary  artefacts. 
Also  it  is  evident  from  the  reconstructed  images  that  there  exist  large  boundary  artefacts  which 
are  due  to  the  incorrect  model.  However,  both  the  quality  and  quantitative  accuracy  of  2D 
image  reconstruction  will  be  improved  by  the  incorporation  of  a  priori  information  as  well  as 
a  multi-step  image  reconstruction  where  regions  of  interest  can  be  identified  and  isolated  for 
regional  reconstruction. 

In  the  case  of  3D  image  reconstruction,  assumptions  regarding  the  breast  not  being 
deformed  might  be  used  and  the  reconstructed  images  provide  relatively  good  absorption 
distributions,  but  the  images  contain  large  artefacts  which  will  likely  confound  their 
interpretation.  The  reason  for  this  can  be  explained  by  considering  the  overall  shape  of 
the  breast  before  and  after  deformation.  The  non-deformed  mesh  at  the  plane  of  measurement 
has  radius  that  is  19  mm  larger  than  the  deformed  mesh.  Additionally,  as  evident  from 
figures  1(a)  and  1(b),  it  is  found  that  the  breast  expands  above,  below  and  between  each  optical 
fibre  once  compressed.  This  large  change  in  shape  of  the  breast  contributes  to  the  large  artefacts 
seen  within  the  reconstruction.  However,  if  one  assumes  that  the  diameter  of  the  breast  is 
known  within  the  measurement  plane,  as  well  as  an  approximate  diameter  above  and  below 
it,  it  is  possible  to  create  a  pseudo-3D  conical  shaped  mesh  for  image  reconstruction.  Images 
reconstructed  using  this  3D  conical  mesh  have  shown  a  better  accuracy  in  recovering  both 
the  absorbing  and  scattering  anomalies.  Finally,  if  accurate  knowledge  regarding  the  breast 
deformation  is  available,  images  of  optical  properties  can  be  reconstructed  which  localize  both 
anomalies  with  great  accuracy,  and  also  contain  little  or  no  artefact  which  otherwise  would 
arise  from  model  mismatch. 

The  goal  of  this  work  is  to  incorporate  this  new  model  of  breast  deformation  with  more 
accurate  information  regarding  the  mechanical  properties  of  the  breast  to  improve  the  NIR 
image  reconstruction.  The  mechanical  property  information  is  readily  available  from  other 
imaging  modalities,  and  the  synthesis  of  this  information  may  provide  fundamentally  new 
information  about  breast  physiologic  response  to  pressure,  and/or  breast  pathology  response 
to  pressure.  An  accurate  model  of  breast  deformation  should  in  principle  allow  us  to  create 
patient  specific  models  and  meshes,  which  would  in-tum  provide  more  clinically  useful  data. 
Work  is  in  progress  to  validate  and  further  improve  the  deformation  model  with  application  to 
the  female  breast  imaging. 
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The  design  and  implementation  of  a  multispectral,  frequency-domain  near  infrared  tomography 
system  is  outlined,  which  operates  in  a  MRI  magnet  for  utilization  of  MR-guided  image 
reconstruction  of  tissue  optical  properties.  Using  long  silica  optical  fiber  bundles,  measurements  of 
light  transmission  through  up  to  12  cm  of  female  breast  tissue  can  be  acquired  simultaneously  with 
MRI  scans.  The  NIR  system  utilizes  six  optical  wavelengths  from  660  to  850  nm  using  intensity 
modulated  diode  lasers  nominally  working  at  100  MHz.  Photomultiplier  tube  detector  gain  levels 
are  electronically  controlled  on  a  time  scale  of  200  ms,  thereby  allowing  rapid  switching  of  the 
source  to  locations  around  the  tissue.  There  are  no  moving  parts  in  the  detection  channels  and  for 
each  source  position,  15  PMTs  operating  in  parallel  allow  sensitivity  down  to  0.5  pW/cm2  at  the 
tissue  surface.  Images  of  breast  tissue  optical  absorption  and  reduced  scattering  coefficients  are 
obtained  using  a  Newton-type  reconstruction  algorithm  to  solve  for  an  optimal  solution  using  the 
measurement  data.  In  medical  imaging,  it  is  beneficial  to  compare  the  same  tissue  volume  as  seen 
by  a  variety  of  modalities,  and  perhaps  more  importantly,  there  is  the  hypothesis  that  one  imaging 
system  which  has  high  spatial  resolution  can  be  used  to  enhance  the  reconstruction  of  another 
system  which  has  good  contrast  resolution.  In  this  study  we  explore  the  synergistic  benefits  of  a 
combined  NIR-MRI  data  set,  specifically  the  ways  in  which  MRI  (i.e.,  high  spatial  resolution) 
enhances  NIR  (i.e.,  high  contrast  resolution)  image  reconstruction.  The  design,  calibration,  and 
performance  of  the  imaging  system  are  described  in  the  context  of  preliminary  phantom  tests  and 
initial  in  vivo  patient  imaging.  Co-registered  MRI  validates  and  improves  optical  property 
estimation  in  2D  tomographic  image  reconstructions  when  specialized  algorithms  are  used.  ©  2004 
American  Institute  of  Physics.  [DOI:  10.1063/1.1819634] 


I.  INTRODUCTION 

Diffuse  optical  tomography  (DOT)  with  near  infrared 
(NIR)  light  can  be  used  to  produce  spatially  resolved  images 
of  tissue  optical  properties.  These  optical  property  maps  can 
be  acquired  at  different  wavelengths  and  combined  to  reveal 
hemoglobin  concentration,  oxygen  saturation,  water  and  fat 
content,  as  well  as  a  description  of  scattering  structures.1-5 
These  latter  parameters  are  important  indicators  of  metabolic 
activity,  functional  processes,  or  presence  and  staging  of  dis¬ 
ease.  NIR  diffuse  tomography  generally  suffers  from  com¬ 
paratively  low  spatial  resolution  due  to  the  multiple  scatter¬ 
ing  events  that  occur  along  each  photon  path.6-8  However, 
the  promise  of  this  imaging  modality  lies  in  the  fact  that  it 
affords  new  physical  bases  for  contrast  in  tissue.  For  ex¬ 
ample,  hemoglobin-based  contrast  in  tumors  relative  to  nor¬ 
mal  tissue  is  exceptionally  high  (i.e,  between  100%-300%).9 
However,  improving  the  limitations  associated  with  its  low 
spatial  resolution  is  fundamental  to  implementing  this  tech¬ 
nology  clinically.  A  priori  knowledge  of  tissue  structure  can 
be  used  to  constrain/guide  the  iterative  NIR  image  recon¬ 
struction  process,  and  improve  the  spatial  resolution  and 
quantitative  accuracy  of  recovered  physiological  parameters. 
Consequently,  NIR  techniques  have  been  combined  with  sev¬ 
eral  high  spatial-resolution,  structure-bearing  imaging  mo¬ 


dalities  including  x-ray  tomosynthesis,10  ultrasound  (US),11 
and  magnetic  resonance  imaging  (MRI),12-14  to  study  human 
tissues  and  small  animals.  In  this  report,  we  present  a  com¬ 
bined  NIR-MRI  system  for  imaging  female  breast  tissue,  and 
explore  the  benefits  of  the  combined  data  set  in  a  planar 
tomographic  geometry  where  the  breast  is  imaged  pendent  in 
a  standard  MR  breast  coil. 

In  recent  years,  several  important  technological  factors 
have  contributed  to  the  advancement  of  NIR  spectral  imag¬ 
ing  applied  to  breast  cancer  characterization,  including  an 
increased  understanding  of  light  and  tissue  interaction,  and 
new  software-based  developments  in  image-reconstruction 
algorithms  for  DOT.5  NIR  radiation  (700-900  nm)  is  nonion¬ 
izing,  light  delivery  and  detection  instrumentation  is  rela¬ 
tively  inexpensive,  and  uncomfortable  tissue  compression  is 
not  required  (as  in  mammographic  level  compression).  De¬ 
tection  systems  with  high  signal-to-noise  ratios  can  routinely 
measure  NIR  light  transmission  through  10  cm  of  breast  tis¬ 
sue.  This  distance  may  increase  to  12  cm  for  breasts  with 
fatty  composition  (i.e.,  lower  radiographic  density  and  NIR 
attenuation).  Photon  propagation  within  the  breast  is  well 
described  by  diffusion  theory  since  the  probability  of  photon 
scattering  is  much  greater  than  absorption.  Light  transmis¬ 
sion  measurements  can  be  combined  with  diffusion  theory  to 
provide  robust  image  reconstruction  of  tissue  optical  coeffi- 
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cients.  We  use  a  Newton-type  algorithm  to  solve  for  the  op¬ 
timal  solution  that  provides  a  minimum  error  between  the 
measured  data  and  predicted  response  from  a  model  of  the 
frequency-domain  diffusion  equation.15  In  theory,  the  ap¬ 
proach  allows  separation  of  the  absorption  and  reduced  scat¬ 
tering  coefficients  (jma  and  /z !s).  Frequency-domain  measure¬ 
ments  of  optical  flux  provide  both  amplitude  and  time-based 
(i.e.,  phase  shift)  information — a  unique  data  set  to  solve  the 
estimation  problem  associated  with  recovering  both  coeffi¬ 
cients  simultaneously. 

A  variety  of  imaging  methods  have  achieved  high- 
spatial-resolution  imaging  with  acceptable  to  excellent  soft 
tissue  contrast,  including  x-ray  computed  tomography  (CT), 
US,  and  MRI.  These  techniques  primarily  provide  images  of 
tissue  structure  and  have  a  limited  ability  to  monitor  param¬ 
eters  related  to  tissue  function  other  than  through  the  intro¬ 
duction  of  exogeneous  contrast  agents.  Alternatively,  nuclear 
medicine  approaches  are  routinely  used  to  image  tissue  func¬ 
tions,  such  as  metabolic  fluorodeoxy glucose  uptake,16  and 
many  commercial  systems  exist  to  co-register  these  images 
with  the  structural  data  derived  from  CT,  US,  and  MRI.  The 
combination  of  high  resolution  structural  imaging  with  lower 
resolution  functional  information  is  a  major  emphasis  in  con¬ 
temporary  medical  imaging,  and  customized  hybrid  imaging 
systems  are  being  developed  to  avoid  the  complications  as¬ 
sociated  with  tissue  movement  between  imaging  exams, 
which  compromises  the  accuracy  of  postprocedure  co¬ 
registration. 

The  application  of  NIR  tomography  to  provide  spatially- 
resolved  functional  information,  such  as  hemoglobin  levels, 
oxygen  saturation,  water,  lipid  and  scatterer  content  will 
likely  be  important,  yet  customized  imaging  systems  which 
couple  to  MRI,  US  or  CT  need  to  be  developed  to  evaluate 
and  exploit  this  potential.  The  pioneering  work  of  Ntziach- 
ristos  et  al.  was  the  first  to  demonstrate  feasibility  of  the 
hybrid  NIR-MR  approach  for  clinical  breast  imaging.  Their 
approach  consisted  of  a  co-planar  array  of  optical  source/ 
sensors  detecting  time-resolved  illuminations  when  applied 
to  the  surface  of  the  breast.  Data  was  presented  as  maps  of 
localized  optical  response  at  each  detector  site  but  was  oth¬ 
erwise  not  processed  into  depth-resolved  images  of  optical 
property  distributions.  Zhu  et  al.n  have  demonstrated  clini¬ 
cal  application  of  a  combined  NIR-US  system.  With  a  hand¬ 
held  probe,  they  simultaneously  measure  reflection  of  ultra¬ 
sound  and  modulated  NIR  excitations  at  the  breast  tissue 
surface,  and  construct  co-registered  images  of  structure  and 
optical  absorption  and  hemoglobin  concentration. 

One  of  the  most  challenging  issues  in  hybrid  system  de¬ 
velopment  is  the  incorporation  of  anatomical  information  as 
prior  constraints  into  NIR  image  reconstruction.  This  has 
been  explored  theoretically  in  several  successful  reports.17-22 
Structural  information  can  be  used  to  guide  functional  image 
reconstruction  through  knowledge  of  tissue  composition  and 
location  by:  (i)  alteration  of  the  objective  functions  used  in 
image  reconstruction,10’23  (ii)  reduction  of  the  number  of  un¬ 
known  parameters  by  treating  regions  of  the  same  tissue  type 
as  single  zones,13,24  or  (iii)  introduction  of  special  regulariza¬ 
tion  schemes  that  can  stabilize  the  inverse  problem  and  em¬ 
phasize  image  contrasts.10,25,26  Optimizing  spatial  resolution 


more  likely  will  depend  on  the  application  of  all  three  tech¬ 
niques,  while  also  being  strongly  influenced  by  the  signal  to 
noise  ratio  of  the  measurements,  the  optical  contrast  avail¬ 
able,  and  the  number  of  projections  used.  Some  investigators 
have  successfully  combined  different  approaches  with  multi- 
stage  image  reconstruction  algorithms.  ’  However,  there  is 
still  no  clear  consensus  on  how  best  to  utilize  the  structural 
information  to  enhance  or  improve  the  recovery  of  functional 
NIR  information. 

This  work  describes  the  first  system  to  combine  multi- 
spectral  frequency-domain  NIR  in  a  planar  tomographic  ge¬ 
ometry  with  MRI  for  imaging  breast  tissue.  NIR  tomography 
has  shown  the  ability  to  localize  changes  in  functional  tissue 
parameters  in  vivo ,  and  MRI  has  the  advantage  of  offering  a 
particularly  rich  amount  of  anatomical  information,  specifi¬ 
cally  about  the  layered  adipose  and  glandular  tissue  structure 
of  the  breast.  This  system  is  designed  to  combine  the  benefits 
of  both  modalities  into  the  construction  of  a  single  quantita¬ 
tive  image.  In  particular,  since  NIR  tomography  has  signifi¬ 
cant  difficulty  with  the  recovery  of  properties  in  layered 
structures,  the  initial  information  from  MRI  can  significantly 
improve  the  estimation  of  breast  properties,  including  in  lo¬ 
calized  regions  such  as  tumors. 

The  design  and  operation  of  the  NIR-MRI  system  ele¬ 
ments  are  described  in  Sec.  II.  The  NIR  component  is  similar 
to  an  imaging  system  described  previously,  which  is  cur¬ 
rently  being  evaluated  clinically,  yet  has  the  unique  design 
feature  of  having  no  moving  parts  in  the  detection  channels, 
allowing  significant  reduction  in  the  NIR  data  acquisition 
time.  In  Sec.  Ill  we  outline  the  theoretical  basis,  and  the 
practical  application  and  utilization  of  image  reconstruction 
in  NIR  tomography.  After  applying  several  approaches  to 
optimizing  this  hybrid  reconstruction,  an  algorithm  is  exam¬ 
ined  that  takes  advantage  of  the  composite  data  set.  In  Sec. 
IV  we  discuss  system  performance  and  present  images  of 
phantom  and  breast  optical  properties  which  are  both  high 
resolution  and  quantitatively  accurate. 

II.  SYSTEM  DESIGN 

This  section  describes  the  four  elements  that  comprise 
the  NIR-MRI  system:  (A)  light  delivery,  (B)  detector  array, 
(C)  fiberoptic  patient  interface,  and  (D)  computer  control  and 
electronics.  Element  (C)  extends  from  (A)  and  (B)  into  a 
MRI  scanner  for  clinical  studies,  as  shown  schematically  in 
Fig.  1.  Photographs  of  the  rack-mounted  system  and  patient 
interface  are  presented  in  Fig.  2.  Section  (E)  briefly  describes 
phantom  fabrication,  and  its  importance  in  imaging  system 
development  and  validation. 

A.  Light  delivery 

The  system  deploys  six  laser  diodes:  661  nm  (40  mW), 
752  nm  (50  mW),  785  nm  (50  mW),  805  nm  (50  mW),  829 
nm  (50  mW),  and  849  nm  (50  mW).  Each  wavelength  is 
amplitude  modulated  at  100  MHz  by  mixing  a  dc  current 
source  (LDX-3220,  ILX  Lightwave,  Bozeman,  MT)  and  an 
ac  current  from  a  frequency  generator  (IFR-2023A,  IFR  Sys¬ 
tems,  Wilmington,  MA),  through  a  bias  T  (#5545,  Picosec¬ 
ond  Pulse  Labs,  Boulder,  CO).  Each  diode  is  held  in  a  laser 


Downloaded  03  Mar  2005  to  129.170.67.67.  Redistribution  subject  to  AIP  license  or  copyright,  see  http://rsi.aip.org/rsi/copyright.jsp 


5264 


Rev.  Sci.  Instrum.,  Vol.  75,  No.  12,  December  2004 


Brooksby  et  al. 


DC  current  source 


6  LOs  on  translation  staff  ^16  WftKa1ed  °P1ical  (13  m) 


AC  (100  MHz) 

®^AC  (100.0005  MHz)U®«— = 


Filter  &  Amp. 


Instrument  Control 
Lock-in  detection 
Image  reconstruction 


FIG.  1.  (Color  online)  Schematic  de¬ 
sign  of  a  dual  modality  NIR-MRI  sys¬ 
tem  (left).  Frequency-domain  NIR  to¬ 
mography  is  performed  inside  the 
MRI  unit  (upper  right).  Six  laser  di¬ 
odes  (660-850  nm)  are  amplitude 
modulated,  and  sixteen  projections 
yield  240  measurements  of  amplitude 
and  phase  of  transmitted  light.  The  fi¬ 
beroptic  array  is  positioned  inside  the 
open  breast  coil,  to  allow  positioning 
along  the  length  of  the  pendant  breast 
(bottom  right). 


tube  (Thorlabs,  Newton,  NJ),  and  mounted  on  a  linear  trans¬ 
lation  stage  (MA2515P5-S2.5  Velmex,  Bloomfield,  NY). 
This  stage  directs  a  specified  wavelength  into  one  of  sixteen 
bifurcated  optical  fiber  bundles  which  were  custom  designed 
for  this  application  (Ceramoptec,  East  Longmeadow,  MA). 
The  248-piece  bundles  (0.37  N.A.,  0.68  packing  fraction)  are 
pure  silica  core  (210  pm),  silicone  clad  (230  pm)  fibers  suit¬ 
able  for  transmission  wavelengths  from  400  nm  to  2400  nm. 
The  source  light  is  delivered  through  the  central  seven  fibers 
in  each  bundle,  and  the  remaining  fibers  surrounding  these 
are  delivered  to  the  detectors.  The  common  end,  which 
makes  contact  with  the  tissue,  has  a  diameter  of  4  mm.  Each 
fiber  bundle  is  13  m  in  length  and  extends  from  the  instru¬ 
ment  cart,  located  outside  of  the  MR  suite,  into  the  bore  of 
the  scanner  (1.5T  whole  body  imager,  GE  Medical  Systems, 
Milwaukee,  WI)  to  the  patient  interface.  The  efficiency  of  the 
optical  switching  is  approximately  50%,  yielding  an  average 
source  power  of  15  mW  at  the  tissue  surface. 

B.  Light  detection 

For  each  source  excitation,  light  transmission  is  recorded 
from  15  surface  locations.  This  signal  is  measured  by  15 


photomultiplier  tubes  (PMT  R6357,  Hamamatsu,  Japan)  op¬ 
erating  in  parallel.  The  gain  of  the  PMTs  is  varied  to  account 
for  the  large  variation  in  light  level  between  detectors 
depending  on  their  distance  from  the  source.  The  gains  are 
set  with  PMT  modules  (HC120,  Hamamatsu)  by  applying 
computer  generated  voltages  between  0.4  and  1.2  V  to  their 
control  lines,  which  sets  the  anode  to  cathode  voltage  be¬ 
tween  approximately  350  V  to  1000  V,  respectively.  Using 
the  higher  gain  settings,  a  PMT  can  reliably  measure  optical 
signals  in  the  pW  range.  The  optimal  gain  levels  are  deter¬ 
mined  prior  to  each  imaging  series.  Each  PMT  is  fixed  to  a 
particular  fiber,  so  it  is  necessary  to  switch  gains  electroni¬ 
cally  during  the  course  of  data  collection.  A  100  M(1  resistor 
was  used  in  the  dynode  chain  of  each  PMT  to  achieve  fast 
settling  times  after  gain  adjustment  (200  ms  for  large  gain 
changes).  Electrical  heterodyning  through  rf  mixers  (Minicir¬ 
cuits,  Brooklyn,  NY)  is  used  to  down  convert  the  100  MHz 
PMT  signal  to  a  lower  frequency  (500  kHz).  This  offset  fre¬ 
quency  is  achieved  with  a  second  frequency- synthesizer 
which  is  synchronized  to  the  one  driving  the  laser  current, 
and  is  set  to  100.0005  MHz.  The  resulting  offset  frequency  is 
filtered  and  amplified  by  a  16  channel  circuit  designed  for 


FIG.  2.  (Color  online)  Photograph  of 
the  rack  mounted,  portable  system 
(left).  System  components  are  marked, 
including  (A)  frequency  generators, 
(B)  optical  switching  stage,  (C)  PMT 
detection  plate,  and  (D)  optical  fibers 
and  patient  interface.  The  fiber-patient 
interface  (at  right)  can  accommodate 
breasts  6-12  cm  in  diameter.  A  close 
up  of  the  tissue  coupling  system  shows 
that  the  fibers  are  spring-loaded  and 
make  light  contact  with  the  full  cir¬ 
cumference  of  a  pendant  breast. 
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this  application  (Audon  Electronics,  Nottingham,  UK),  then 
read  by  the  computer.  Lock-in  detection  is  executed  in  soft¬ 
ware  to  extract  amplitude  and  phase  data  for  each  of  the 
detectors  in  parallel. 


C.  Fiberoptic  patient  interface 

The  MR  exam  is  performed  using  a  breast  array  coil 
(MRI  Devices,  Waukesha,  WI)  that  offers  high-resolution 
imaging.  The  coil  also  provides  an  open  architecture,  which 
allows  for  the  integration  of  the  NIR-breast  interface  shown 
in  Fig.  1.  A  circular  ring  machined  from  polyvinyl  chloride 
(PVC)  positions  the  common  ends  of  the  sixteen  fibers 
around  the  full  circumference  of  the  pendant  breast.  Phos¬ 
phor  bronze  compression  springs  (k=0.09  lbs/in.,  Ace  Wire 
Spring  and  Form,  McKees  Rocks,  PA)  guide  each  fiber 
through  holes  in  the  ring,  into  light  contact  with  the  tissue 
surface.  The  ring  separates  into  equal  halves  so  that  it  can 
easily  be  moved  from  one  breast  to  the  other.  The  angular 
separation  between  each  fiber  is  21°,  except  between  fibers 
adjacent  to  the  line  of  ring  de vision,  where  the  separation  is 
33°.  The  ring  can  be  positioned  vertically,  such  that  the  plane 
of  measurements  intersects  the  region  of  interest  in  tissue. 
This  design  allows  each  fiber  to  move  independently,  there¬ 
fore  their  radial  positions  may  vary  by  several  millimeters.  In 
order  to  avoid  serious  artifacts  in  reconstructed  images,  the 
position  of  each  fiber  must  be  accurately  known  during  data 
acquisition.  Annular  fiducial  markers  (MM  3005,  IZI  Medi¬ 
cal  Products,  Baltimore,  MD)  are  fixed  to  each  fiber  and  can 
be  located  with  submillimeter  accuracy  in  the  MRI. 


D.  Computer  system 

A  PC  running  Lab  view  software  (National  Instruments, 
Austin,  TX)  is  used  to  control  all  light  delivery  and  detection 
equipment.  The  laser  current  source  and  frequency  generator 
parameters  are  set  by  a  general  purpose  interface  bus  (GPIB, 
NI).  The  linear  translation  stage  is  addressed  through  the 
serial  port.  An  analog  output  board  (NI)  is  used  for  PMT  gain 
control.  A  multipurpose  data  acquisition  (DAQ)  board  ac¬ 
quires  the  16  analog  input  channels  and  the  single  reference 
channel.  This  board  also  provides  six  digital  output  lines  to 
the  high  power  radio-frequency  switch  for  the  laser  sources. 
For  each  source  position,  15  signals  from  the  detector  system 
are  amplified  by  a  gain  of  1000  and  low  pass  filtered  to 
prevent  aliasing  prior  to  the  DAQ  board  using  a  16  channel 
amplifier  and  filter  network  mounted  in  a  BNC  coupled  box 
(Audon  Electronics,  Nottingham,  UK).  Data  are  acquired  for 
500  ms,  and  phase  and  amplitude  of  each  signal  are  calcu¬ 
lated  and  written  to  file.  Including  the  time  required  to  deter¬ 
mine  optimal  gain  values,  measurement  with  6  wavelengths 
takes  approximately  4  min.  The  MR  exam  is  controlled  sepa¬ 
rately,  operated  in  parallel,  and  a  full  volume  breast  MRI  is 
of  similar  duration.  A  FORTRAN,  or  MATLAB  based  recon- 
struction  program  reads  and  processes  the  NIR  data.  MR 
images  are  processed  off-line  with  an  addition  to  the  MATLAB 
software  package,  and  incorporated  into  a  modified  iterative 
optical  property  reconstruction  (see  Sec.  III). 


E.  Phantom  design 

Tissue- simulating  phantoms  with  known  property  distri¬ 
butions,  geometries,  and  imaging  orientations  are  commonly 
used  to  validate  imaging  systems.  We  have  developed  a 
recipe  for  producing  gelatin  materials  with  desired  optical 
properties,  and  a  shelf  life  of  several  months.31  A  heated 
mixture  of  water,  gelatin  (G2625,  Sigma  Inc.),  India  ink  (for 
absorption),  and  titanium  dioxide  powder  (for  scatter)  (Ti02, 
Sigma  Inc.)  is  poured  into  a  mold  of  a  desired  shape,  and 
solidified  by  cooling  to  room  temperature.  Variation  in  the 
water  concentration  provides  MR  contrast,  and  variable  gel 
stiffness.  The  phantom  imaged  here  (Sec.  IV  B)  combines 
three  gels  with  different  optical  properties  in  an  irregular 
structure. 

III.  DATA  PROCESSING  AND  IMAGE 
RECONSTRUCTION 

Quantitative  NIR  imaging  with  model-based  methods  re¬ 
quires  (A)  important  instrument  calibration  procedures,  and 
(B)  a  reconstruction  algorithm  that  incorporates  an  accurate 
model  of  light  propagation  in  tissue. 

A.  System  calibration 

Calibration  issues  and  other  practical  considerations  as¬ 
sociated  with  our  NIR  imaging  approach  have  been  dis- 
cussed  in  detail  elsewhere.  Two  important  procedures  are 
briefly  noted  here:  (1)  detector  calibration,  and  (2)  homoge¬ 
neous  phantom  calibration.  First,  the  amplitude  and  phase 
response  of  each  detection  channel  must  be  characterized  in 
order  to  remove  systematic  noise  in  the  data  acquisition 
hardware.  Each  detector  is  exposed  to  the  same  optical  sig¬ 
nal,  and  the  differences  in  log  amplitude  and  phase  are  used 
as  correction  factors.  The  log  amplitude  response  of  the  PMT 
is  plotted  against  the  log  of  the  input  power  for  each  gain 
setting.  A  log-log  regression  is  performed  and  the  coeffi¬ 
cients  are  used  to  calibrate  detected  PMT  amplitude  in  terms 
of  optical  power.  The  phase  does  not  fluctuate  significantly 
with  changing  light  level  for  a  single  gain  setting  (i.e.,  mini¬ 
mal  phase-amplitude  cross-talk),  but  is  altered  dramatically 
with  changing  gain.  Relative  phase  differences  between  de¬ 
tectors  are  stored  for  calibration.  These  calibration  curves  are 
very  similar  to  those  created  by  McBride  et  al.21  This  char¬ 
acterization  needs  to  be  performed  only  once  as  long  as  the 
system  is  not  modified. 

The  second  important  practical  procedure  is  the  correc¬ 
tion  for  inter  fiber  variations  and  coupling  issues,  which  is 
accomplished  through  a  homogeneous  phantom  calibration 
process.  ’  This  accounts  for  offsets  due  to  optical  fiber 
differences  in  transmission  and  alignment,  as  well  as  any 
errors  in  discretization  or  data-model  mismatch.  A  homoge¬ 
neous  phantom  is  generally  measured  each  day,  and  after 
system  changes.  The  differences  between  data  measured 
from  the  phantom,  and  data  calculated  from  the  model  are 
stored  and  subtracted  from  measurements  of  the  heteroge¬ 
neous  phantom  or  tissue  under  investigation.  A  homogeneous 
fitting  algorithm  is  used  to  determine  the  and  jj!s  values 
supplied  to  the  model  calculation.  This  algorithm  can  also  be 
used  to  calculate  the  initial  optical  properties  specified  in 
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iterative  reconstruction  of  heterogeneous  media.  When  deal¬ 
ing  with  tissues  having  arbitrary  shape,  the  effectiveness  of 
this  fitting  algorithm  and  homogeneous  phantom  calibration 
hinges  on  the  accurate  specification  of  source  and  detector 
locations.  The  ability  to  extract  accurate  fiber  positions  from 
MRI  scans  preserves  the  integrity  of  this  method  for  nonuni¬ 
form  boundary  data. 

B.  FEM  analysis 

Data  acquired  from  the  detection  system  is  processed  by 
a  FEM  based  reconstruction  algorithm  to  generate  tomogra¬ 
phic  images  of  absorption  and  reduced  scattering  coefficients 
simultaneously.  The  algorithm  exploits  the  frequency- 

domain  diffusion  equation  approximation  to  light  behavior  in 

28 

a  highly  scattering  medium, 

(  too  \ 

^a(r)  +  “  J<I>(r,a>)  =  S(r,(o), 

(1) 

where  S(r,(o)  is  an  isotropic  light  source  at  position  r, 
<F(r,w)  is  photon  density  at  r,  c  is  the  speed  of  light  in 
tissue,  o)  is  the  frequency  of  light  modulation,  [ia  is  the  ab¬ 
sorption  coefficient,  and  D=  l/[3(/zfl+/^)]  is  the  diffusion 
coefficient.  The  reduced  (transport)  scattering  coefficient  is 
given  by  jul's  =  /x^(l-g),  where  g  is  the  mean  cosine  of  the 
single  scatter  function  (the  anisotropy  factor),  and  juls  is  the 
scattering  coefficient. 

For  a  given  fxa  and  /uls  distribution,  the  diffusion  equation 
is  used  to  predict  the  optical  flux  at  the  detector  sites  for  each 
source  excitation.  In  the  inverse  problem  (image  reconstruc¬ 
tion),  the  goal  is  the  estimation  of  optical  properties  at  each 
FEM  node,  based  on  measurements  of  optical  flux  at  the 
detector  sites  on  the  tissue  surface.  This  is  achieved  numeri¬ 
cally  by  minimizing  the  difference  between  the  calculated 
data  Oc,  and  measured  data,  Om,  for  all  source/detector 
combinations  (NM).  Typically, 

NM 

V=2(3>f-4f)2  (2) 

i= 1 

is  minimized  in  a  least  squares  sense  by  setting  the  first  de¬ 
rivative  equal  to  zero,  and  using  a  Newton-Raphson  ap¬ 
proach  to  find  the  set  of  optical  property  values  which  ap¬ 
proximate  the  point  of  stationarity.  We  use  a  Levenberg 
Marquardt  algorithm,  and  repeatedly  solve  the  equation 

a  =  (JTJ+XI)~1JTb ,  (3) 

where  Z?  =  (Oc-T>M)r  is  our  data  vector,  and  a  is  the  solution 
update  vector,  a-[SDj\  Sjmajl  defining  the  difference  between 
the  true  and  estimated  optical  properties  at  each  recon¬ 
structed  nodey.  Here,  X  is  a  regularization  factor  to  stabilize 
matrix  inversion  and  J  is  the  Jacobian  matrix  for  our  model, 
which  is  calculated  using  the  Adjoint  method.33 

Improving  NIR  reconstructions  by  incorporating  MRI 
data  has  been  explored  in  previous  work,  ’  and  by  other 
authors.  ’  ’  Techniques  used  to  produce  the  images  shown 
in  Sec.  IV  include:  (i)  accurately  defining  the  imaging  vol¬ 
ume,  (ii)  tailoring  a  regularization  scheme  which  optimizes 
the  reconstructed  contrast  of  a  suspicious  area  in  the  image, 
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FIG.  3.  (Color  online)  Repeatability  assessment  for  (a)  log  amplitude  and 
(b)  phase  of  NIR  component  of  combined  NIR-MRI  system  and  stand-alone 
NIR  tomography  system  at  Dartmouth.  The  performance  of  the  two  systems 
is  comparable.  In  routine  operation,  PMT  voltage  signals  are  above  1.0 

X  1(T5  V. 

and  (iii)  reducing  the  number  of  unknown  parameters  by 
segmenting  tissue  types  visible  in  the  MRI.  Defining  the  im¬ 
aging  volume,  (i),  is  relatively  straight  forward  and  can  be 
accomplished  by  creating  a  structured  finite  element  mesh 
from  the  MRI.  This  mesh  will  often  contain  regional  differ¬ 
ences  depending  on  tissue  types  present,  and  an  irregular 
outer  boundary  due  to  impressions  caused  by  fiber-tissue 
contact.  In  simulation  studies,  not  presented  here,  it  has  been 
found  that  image  reconstruction  accuracy  is  easily  degraded 
if  the  mesh  (2D  or  3D)  does  not  represent  the  true  outer 
tissue  boundary.  Step  (ii)  can  be  accomplished  by  changing 

XI  in  Eq.  (3)  to  XA ,  where  A  is  a  regularization  matrix,  or 
filter  matrix.  Regularization  can  be  thought  of  as  a  smooth¬ 
ing  operator,  where  one  can  apply  selective  smoothing  by 
linking  together  the  property  updates  for  all  nodes  associated 
with  the  same  region  or  tissue  type.  Modifications  to  the 
Jacobian  matrix  size  in  a  parameter  reduction  technique  are 
used  to  implement  step  (iii).23 


IV.  PERFORMANCE  AND  EXPERIMENTAL  RESULTS 

In  this  section  we  (A)  compare  measurements  from  this 
system  to  those  from  a  similar  NIR  tomography  instrument 
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TABLE  I.  Optical  properties  at  785  nm  for  15  tissue  simulating  phantoms, 
and  %  differences  from  values  determined  using  another  NIR  tomography 
system. 


/A( mm  9 
(NIR-MRI) 

%  difference 

rfmm  9 

(NIR-MRI) 

%  difference 

Mean 

0.0055 

-6.2 

1.34 

9.1 

Max 

0.0102 

-23.4 

1.91 

29.6 

Std.  Dev. 

0.0024 

8.7 

0.41 

7.9 

described  previously27  for  a  number  of  tissue  simulating 
phantoms.  NIR-MRI  phantom  studies  are  described  in  (B), 
and  in  vivo  images  are  shown  in  (C).  For  all  results  presented 
here,  we  used  a  single  laser  diode  (785  nm)  for  convenience, 
and  2D  modeling  and  image  reconstruction.  System  perfor¬ 
mance  and  image  quality  at  the  other  five  wavelengths  are 
comparable,  and  3D  imaging  is  readily  achievable. 

A.  System  performance 

Measurement  repeatability  in  terms  of  phase  and  log  am¬ 
plitude  error  was  assessed  by  serially  imaging  a  phantom 
with  optical  properties  similar  to  those  of  average  breast  tis¬ 
sue  (juia= 0.004,  /z^=1.35).  The  average  rms  error  at  each 
detector  site  was  determined  to  be  0.26%  in  ac  intensity  and 
1.04°  in  phase.  These  values  compare  with  those  obtained 
from  the  system  described  by  McBride  et  al.  (0.32%  in  ac 
intensity  and  0.48°  in  phase),  rms  error  for  each  of  the  240 
source-detector  combinations  for  both  systems  is  plotted  ver¬ 
sus  the  PMT  signal  in  Fig.  3.  For  both  phase  and  amplitude, 
error  sharply  increases  when  incident  light  falls  below  ap¬ 
proximately  0.5  pW.  These  points  are  excluded  in  the  calcu¬ 
lation  of  average  rms  error.  As  an  additional  comparison,  we 
used  both  devices  to  measure  a  collection  of  homogeneous 
phantoms  (A=15)  of  varying  diameters  (73-91  mm) 
and  optical  properties  (/z^  =  0.0023-0.0 102  mm-1,  /j!s 
=  0.33-1.91  mm-1).  After  processing  the  measurements  with 
the  calibration  procedure  described  in  Sec.  Ill,  we  used  the 
homogeneous  fitting  algorithm  to  determine  a  single  /Jia  and 
ju!s  for  each  material.  Table  I  shows  the  optical  coefficients 
obtained  with  the  NIR-MRI  system,  along  with  a  measure  of 
their  discrepancy  with  those  obtained  with  our  stand-alone 
DOT  instrumentation.  We  observed  good  agreement  between 


TABLE  II.  Approximate  optical  properties  of  gelatin  phantom. 


Outer  layer 

Inner  layer 

Inclusion 

ft(mm  9 

0.0044 

0.0062 

0.02 

^(mm-9 

0.51 

0.68 

0.9 

the  two  systems,  consistent  throughout  the  range  of  phantom 
properties.  The  absorption  and  scattering  coefficients 
show  correlation  coefficients  of  R2  =  0.984  and  R2= 0.980, 
respectively. 

B.  Phantom  imaging 

Due  to  the  limited  spatial  resolution  of  DOT,  layered 
media,  small  objects,  and  low  contrast  heterogeneity  pose 
key  challenges  in  image  reconstruction.  The  capability  of  the 
presented  system  to  address  these  challenges  was  investi¬ 
gated  by  imaging  a  phantom  comprised  of  three  gels  with 
different  optical  properties.  The  phantom  is  cylindrical  and 
the  boundary  between  the  outer  layer  and  inner  layer  is  ir¬ 
regular.  A  two  centimeter  spherical  inclusion  is  embedded 
within  the  inner  layer.  The  optical  coefficients  of  each  gel  are 
known  (Table  II),  as  each  material  was  created  using  a  prac¬ 
ticed  recipe  to  give  desired  values.  Furthermore,  the  true 
value  of  each  layer  was  validated  by  creating  a  separate  ho¬ 
mogenous  phantom  for  each  layer  (at  the  same  time  as  cre¬ 
ating  layered  phantom)  and  measuring  bulk  properties  with  a 
homogeneous  fitting  algorithm.  To  increase  MRI  contrast  be¬ 
tween  layers,  Omniscan™  (gadodiamide)  was  added  to  the 
inner  layer  (0.005  g/ ml). 

A  photograph  of  the  phantom,  alongside  another  spheri¬ 
cal  inclusion  is  shown  in  Fig.  4(a).  Figure  4(b)  shows  a  Tl- 
weighted,  gradient  echo  MRI  (25  ms  TR,  3  ms  TE,  45°  flip 
angle)  crosssection  of  the  phantom  in  the  plane  of  the  optical 
fibers.  This  was  used  to  create  a  2D  structured  finite  element 
mesh  and  to  locate  the  positions  of  the  16  fibers  [Fig.  4(c)]. 
Fiber  fiducials  were  not  used  in  this  experiment,  but  the  16 
impressions  caused  by  each  optical  fiber  are  clearly  visible 
around  the  perimeter  of  the  gel.  Each  of  the  three  regions  are 
also  visible,  corresponding  to  each  of  the  three  types  of  op¬ 
tically  variant  gel.  3D  meshing  and  imaging  is  also  possible, 
given  that  a  stack  of  MR  slices  represent  the  full  volume 


FIG.  4.  (Color  online)  (a)  The  gelatin  phantom,  (b)  T1 -weighted  MRI  image  of  the  structure  of  the  phantom  cross  section.  The  inner  layer  (light  region) 
contains  gadodiamide  for  MR  contrast  (0.005  g/ml).  The  inner  inclusion  is  a  2  cm  diameter  sphere,  (c)  Finite  element  mesh  derived  from  the  three-layered 
structure  image  in  (b).  The  circles  on  the  outer  boundary  indicate  the  fiber  locations.  Axes  are  in  millimeters. 
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FIG.  5.  (Color  online)  Images  in  (a) 
show  absorption  coefficient  (left)  and 
reduced  scattering  coefficient  (right) 
reconstructions  without  MRI  guid¬ 
ance.  In  (b),  reconstructions  improve 
when  the  interior  structural  informa¬ 
tion  of  the  MRI  is  incorporated.  Com¬ 
paring  these  two  images,  the  latter  has 
a  reduction  in  artifacts  in  the  reduced 
scattering  coefficient  image,  and  both 
the  spatial  resolution  and  contrast  have 
improved.  The  recovered  values  com¬ 
pare  favorably  with  the  approximate 
true  values  shown  in  Table  II. 


structure  of  the  phantom.  The  system  design  allows  for  NIR 
data  acquisition  in  multiple  planes  if  desired. 

Using  log  amplitude  and  phase  data,  images  of  optical 
properties  were  reconstructed  with  two  different  algorithms. 
The  first  algorithm  solves  Eq.  (3)  without  a  priori  guidance, 
except  for  the  use  of  the  mesh  and  source  locations  from  Fig. 
4(c).  The  corresponding  reconstructed  images  are  presented 
in  Fig.  5(a).  The  second  algorithm  uses  the  layered  structure 
from  the  MRI  to  constrain  Eq.  (3).  The  regularization  param¬ 
eter  associated  with  each  reconstructed  node  is  adjusted 
based  upon  its  location  and  area  of  influence.  Eower  values 
are  used  for  small  regions  close  to  the  center  of  the  recon¬ 
structed  volume,  whereas  peripheral  regions  (which  are 
prone  to  large  artifacts)  have  larger  regularization.  A  regular¬ 
ization  matrix  (A),  or  filter  matrix,  based  upon  the  a  priori 
structural  information  from  the  MRI  was  also  used  to  further 
improve  the  algorithm.  The  optimal  values  for  X  and  A  were 
determined  in  simulation  studies  of  this  geometry,  using 
similar  contrast  and  noise  levels.  These  can  also  be  chosen 
automatically  once  an  empirical  knowledge  of  their  effect  is 
established.  The  images  reconstructed  from  this  modified 
constrained  algorithm  are  shown  in  Fig.  5(b). 

Qualitatively,  the  algorithm  that  uses  MRI  information  to 
guide  the  iterative  process  performs  much  better.  The  optical 
property  images  in  Fig.  5(a)  are  blurry  and  “edge  artifacts” 
are  clearly  visible,  especially  for  p!s.  The  absorbing  sphere 
near  the  center  of  the  phantom  is  recovered  with  poor  spatial 
resolution,  and  its  contrast  is  underestimated.  The  property 
maps  in  Fig.  5(b)  have  improved  resolution.  Spatial  resolu¬ 
tion  and  contrast  of  the  spherical  inclusion  are  better.  As  a 
quantitative  measure  of  the  accuracy  of  image  reconstruc¬ 
tion,  we  compute  the  rms  error  between  the  target  image  and 
the  reconstructions,  node  by  node.  This  measure  indicates  a 
dramatic  improvement  for  the  second  method  [rms(/zj 


=  0.0024,  rms(/^)  =  0.19]  relative  to  the  first  [rms(/zj 
=  0.0034,  rms(/z')  =  0.25]. 

C.  In  vivo  imaging 

The  Institutional  Review  Board  (IRB)  at  the  Dartmouth 
Hitchcock  Medical  Center  approved  this  clinical  examina¬ 
tion  protocol,  and  written  informed  consent  was  obtained 
from  all  volunteering  women.  Figures  6(a)  and  6(b)  show  a 
photograph  of  a  subject  lying  on  the  examination  platform, 
and  an  anatomical  axial  MR  image  through  the  breast.  NIR 
measurements  (4  wavelengths  in  this  case)  and  a  full  volume 
MRI  (50  coronal  slices,  25  ms  TR,  6  ms  TE,  45°  flip  angle, 
2  mm  slice  thickness)  were  obtained  in  less  than  10  min.  Fig. 
7(a)  shows  an  anatomically  coronal  MRI  with  16  fiducial- 
marked  fibers  (appearing  as  bright  white  spots  outside  the 
breast).  The  radiographic  density  of  this  participant  is  hetero¬ 
geneously  dense  (HD),  and  a  large  interior  region  of  glandu¬ 
lar  tissue  is  easily  defined  in  the  FEM  mesh  shown  in  Fig. 
7(b). 

As  with  the  phantom  study,  we  present  images  of  p,a  and 
p!s  at  785  nm  reconstructed  with  two  different  algorithms. 
The  first  result,  obtained  by  solving  Eq.  (3),  without  a  priori 
guidance,  is  shown  in  Fig.  7(c).  As  expected,  we  see  higher 
absorption  and  scatter  in  the  glandular  region  (central)  rela¬ 
tive  to  the  adipose  tissue  (peripheral).  However,  the  region  of 
increase  does  not  span  the  full  area  expected,  and  heteroge¬ 
neity  is  visible  (especially  around  the  perimeter  of  the  im¬ 
age).  The  second  algorithm  assumes  homogeneous  optical 
properties  for  each  tissue  type,  and  utilizes  parameter 
reduction,  which  leads  to  a  “fitting”  for  four  values: 

AU  adipose- 0-003,  fls  adipose-  0-93,  AU  glandular- 0-006, 

pls  glandular"  1-12  [Fig.  7(d)].  This  algorithm  is  robust  to  noise 
and  converges  after  a  few  iterations.  The  result  is  quantita- 


FIG.  6.  (Color  online)  Photograph  in 
(a)  shows  a  female  volunteer  prepared 
for  the  simultaneous  MRI-NIR  exam. 
In  an  anatomically  axial  T1 -weighted 
MRI  slice  from  the  right  breast,  (b), 
fiducial  markers  (outside  the  breast) 
indicate  the  location  of  the  fiber  plane. 
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FIG.  7.  (Color  online)  In  (a),  an  anatomically  coronal  Tl-weighted  MRI  displays  adipose  (outer)  and  glandular  (inner)  tissue  types.  A  two  layer  structured 
FEM  mesh  and  source  locations,  created  from  (a),  is  shown  in  (b).  Absorption  and  reduced  scattering  coefficient  (mm-1)  reconstructions  in  (c)  are  obtained 
without  utilizing  the  internal  structure  of  (b).  In  (d),  MRI  data  guides  a  two-region  parameter  fitting  algorithm.  Relative  to  (c),  resolution  has  improved  and 
contrast  has  increased — showing  higher  absorption  and  scatter  in  glandular  relative  to  adipose  tissue. 


tively  logical,  and  similar  values  are  recovered  using  spa¬ 
tially  varying  regularization  (as  described  in  the  phantom 
reconstruction). 

V.  DISCUSSION 

This  article  describes  a  novel  in  vivo  breast  imaging  sys¬ 
tem  that  synergistically  combines  NIR  tomography  with 
MRI.  Tissue  structures  visible  with  high  resolution  in  MRI 
can  be  applied  a  priori  to  optical  property  reconstructions 
from  frequency-domain  NIR  measurements.  Thus,  the  recon¬ 
struction  process  can  be  optimized  to  produce  high  resolu¬ 
tion,  quantitatively  accurate  maps  of  absorption  and  reduced 
scattering  coefficients,  and  ultimately  physiologically  rel¬ 
evant  parameters.  Various  physiological  tissue  types  exhibit 
significant  contrast  in  the  NIR,  primarily  though,  the  com¬ 
bined  system  could  be  very  effective  at  locating  and  diagnos¬ 
ing  breast  tumors. 

The  NIR  component  provides  multispectral  (6  wave¬ 
lengths)  frequency-domain  (log  amplitude  and  phase)  data 
from  16  fiber-tissue  contact  positions  around  the  breast’s  cir¬ 
cumference.  Data  acquisition  is  fully  automated,  and  a  com¬ 
plete  set  of  measurements  (240  source-detector  pairs)  for  one 
wavelength  requires  approximately  40  s.  Light  detection  is 
highly  sensitive  (subpicowatt  limit)  with  low  noise  (rms 
<0.26%,  <1.04°  in  phase).  All  optical  elements  and  controls 
are  mounted  in  a  portable  cart,  and  operation  inside  strong 


magnetic  fields  is  facilitated  with  long  optical  fibers  and  an 
easy-to-use  positioning  system/patient  interface. 

The  deployment  of  dual  modality  NIR  imaging  systems 
in  clinical  applications  has  been  limited  to  date,  mainly  be¬ 
cause  of  the  complexity  of  the  image  reconstruction  problem. 
Here,  a  representative  phantom  and  in  vivo  study  using  one 
wavelength  (785  nm)  are  presented.  We  have  shown  that 
co-registered  MRI  validates  and  improves  optical  property 
estimation  in  2D  tomographic  image  reconstructions  when 
specialized  algorithms  are  used.  Future  work  will  involve  3D 
modeling  and  reconstruction,  which  could  further  improve 
both  qualitative  and  quantitative  aspects  of  the  recovered  co¬ 
efficient  values. 

Preliminary  results  are  encouraging,  and  have  allowed  us 
to  optimize  reconstruction  techniques,  and  automate  con¬ 
straint  selection.  We  have  performed  several  phantom  stud¬ 
ies,  and  demonstrated  the  feasibility  of  imaging  volunteers 
with  healthy  breasts.  Through  the  study  of  more  healthy  sub¬ 
jects,  with  different  radiographic  densities,  we  aim  to  com¬ 
pare  functional  parameters  of  adipose  versus  glandular  tis¬ 
sue.  To  test  the  system’s  ability  at  diagnosing  tumors,  cancer 
patients  will  be  recruited,  and  the  simultaneous  exam  will 
likely  involve  MR  contrast  enhancement. 
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Effects  of  refractive  index  on  near-infrared 
tomography  of  the  breast 


Hamid  Dehghani,  Ben  A.  Brooksby,  Brian  W.  Pogue,  and  Keith  D.  Paulsen 


Near  infrared  (NIR)  optical  tomography  is  an  imaging  technique  in  which  internal  images  of  optical 
properties  are  reconstructed  with  the  boundary  measurements  of  light  propagation  through  the  medium. 
Recent  advances  in  instrumentation  and  theory  have  led  to  the  use  of  this  method  for  the  detection  and 
characterization  of  tumors  within  the  female  breast  tissue.  Most  image  reconstruction  approaches  have 
used  the  diffusion  approximation  and  have  assumed  that  the  refractive  index  of  the  breast  is  constant, 
with  a  bulk  value  of  approximately  1.4.  We  have  applied  a  previously  reported  modified  diffusion 
approximation,  in  which  the  refractive  index  for  different  tissues  can  be  modeled.  The  model  was  used  to 
generate  NIR  data  from  a  realistic  breast  geometry  containing  a  localized  anomaly.  Using  this  simulated 
data,  we  have  reconstructed  optical  images,  both  with  and  without  correct  knowledge  of  the  refractive- 
index  distribution  to  show  that  the  modified  diffusion  approximation  can  accurately  recover  the  anomaly 
given  a  priori  knowledge  of  refractive  index.  But  using  a  reconstruction  algorithm  without  the  use  of 
correct  a  priori  information  regarding  the  refractive-index  distribution  is  shown  as  recovering  the 
anomaly  but  with  a  degraded  quality,  depending  on  the  degree  of  refractive  index  mismatch.  The  results 
suggest  that  provided  the  refractive  index  of  breast  tissue  is  approximately  1.3-1. 4,  their  exclusion  will 
have  minimal  effect  on  the  reconstructed  images.  ©  2005  Optical  Society  of  America 
OCIS  codes:  170.3010,  170.3830,  170.6960. 


1.  Introduction 

Near-infrared  (NIR)  optical  tomography,  is  a  nonin- 
vasive  imaging  modality  in  which  images  of  internal 
optical  properties  are  reconstructed  by  use  of  mea¬ 
sured  transmission  data.1-8  More  specifically,  a  num¬ 
ber  of  optical  fibers  are  placed  around  the  surface  of 
the  tissue  to  be  imaged  and  light  is  transmitted  from 
one  fiber  at  a  time,  while  all  other  fibers  are  used  to 
measure  the  exiting  light  at  discrete  locations  around 
the  tissue.  The  recovery  of  the  interior  optical  prop¬ 
erties  is  treated  as  a  nonlinear  inverse  problem  based 
on  measurements  at  the  exterior  of  the  domain  to  be 
imaged.  The  light  source  can  be  operated  in 
continuous-wave  (cw)  pulsed  or  intensity-modulated 
modes,  giving  rise  to  output  measurements  of  light 
intensity,  temporal  point-spread  function,  or  phase 
and  amplitude,  depending  on  the  strategy  adopted. 


The  authors  are  with  the  Thayer  School  of  Engineering,  Dart¬ 
mouth  College,  8000  Cummings  Hall,  Hanover,  New  Hampshire 
03755-8000.  H.  Dehghani’s  e-mail  address  is  hamid.dehghani@ 
dartmouth.edu. 

Received  26  July  2004;  revised  manuscript  received  9  November 
2004;  accepted  17  November  2004. 

0003-6935/05/101870-09$15.00/0 

©  2005  Optical  Society  of  America 


Typically  light  sources  are  applied  at  a  number  of 
NIR  wavelengths,  each,  in  turn,  giving  rise  to  addi¬ 
tional  spectrally  independent  boundary  measure¬ 
ments.  The  time-resolved  or  frequency-domain  data 
allows  for  reconstruction  of  internal  optical  property 
distributions  of  absorption  (p,a)  and/or  reduced  scat¬ 
ter  (p,s')  for  each  wavelength,  which  can  in  turn  be 
used  to  calculate  images  of  the  concentrations  of  the 
dominant  chromophores,  including  oxyhemoglobin, 
deoxyhemoglobin  and  water,  as  well  as  scattering 
power  and  amplitude.9  These  images  are  the  end  re¬ 
sults  of  NIR  tomography,  and  they  may  complement 
other  imaging  modalities,  potentially  defining  NIR 
tomography  as  a  functional  imaging  technique. 

For  image  reconstruction  with  the  measured 
boundary  data,  several  different  algorithms  can  be 
used,  depending  on  datatype,  number  of  measure¬ 
ments,  and  domain  under  investigation.10  Many  im¬ 
aging  algorithms  rely  on  a  numerical  model,  whereby 
the  process  for  iteratively  updating  the  internal  dis¬ 
tribution  of  absorption  and  scatter  is  to  match  the 
measured  data  to  model-based  predictions.  Assuming 
that  the  tissue  scatter  dominates  absorption  and  that 
the  imaging  domain  is  large  enough  so  that  the  small¬ 
est  source- detector  fiber  distance  is  greater  than  a 
few  mean  scattering  lengths,  the  diffusion  approxi¬ 
mation  (DA)  models  the  light  propagation  effectively. 


1870 


APPLIED  OPTICS  /  Vol.  44,  No.  10  /  1  April  2005 


Recently,  advances  have  been  made  in  represent¬ 
ing  bulk  internal  refractive-index  (RI)  variation  in 
tissue  within  the  numerical  model.11-13  Jiang  and 
Xu12  have  used  a  higher-order  DA  that  also  takes  into 
account  the  RI  as  a  spatially  varying  property  and 
have  presented  reconstructed  images  of  internal  ab¬ 
sorption,  reduced  scatter,  and  refractive  index  using 
cw  data.  Our  earlier  study  showed  the  implementa¬ 
tion  of  RI  spatial  variation  within  the  DA  using  a 
finite-element  modeling  approach.11  Specifically,  we 
defined  an  internal  boundary  condition  that  allows 
for  discontinuous  internal  variations  of  RI  between 
regions  of  contrast  (or  slowly  varying)  optical  proper¬ 
ties.  This  approach  provided  a  good  match  between 
three-dimensional  finite-element  model  results  and 
Monte  Carlo  simulations,  as  well  as  controlled  exper¬ 
imental  measurements,14  suggesting  that  the  tech¬ 
nique  is  a  valid  approximation. 

The  absolute  values  of  RI  for  breast  tissue  types  are 
still  a  subject  under  investigation.  Estimates  are  dif¬ 
ficult  to  obtain  accurately  in  vivo ,  but  some  data  have 
been  reported  for  various  tissue  types.14 15  Although 
adipose  tissue  (fatty  layer)  has  been  measured  to  be 
1.455,  to  the  best  of  our  knowledge  no  results  exist  for 
fibroglandular  tissue.  Nonetheless,  the  RI  of  glandu¬ 
lar  tissue  is  believed  to  be  lower  than  that  of  adipose 
(e.g.,  values  of  1.4  have  been  assumed15’16).  The  ra¬ 
tionale  for  a  lower  value  is  sound  in  that  estimates  of 
the  composition  of  fibroglandular  tissue  place  its  wa¬ 
ter  (typically  greater  than  60%)  and  blood  content 
(greater  than  1%)  to  be  high,  indicating  that  its  bulk 
RI  is  likely  close  to  that  of  water.17 18 

In  this  paper,  we  have  investigated  at  the  effect  of 
variation  of  RI  within  a  realistic  two-dimensional 
model  of  a  female  breast,  by  assuming  distinct  values 
for  the  adipose  and  glandular  tissue  compartments. 
Using  our  current  dual  modality  magnetic  resonance 
imaging  (MRI)-NIR  system,  we  can  collect  NIR  and 
MRI  data  simultaneously19  allowing  accurate  a  priori 
information  regarding  the  tissue  structure  under  in¬ 
vestigation.  Specifically,  we  have  generated  NIR  data 
from  a  two-dimensional  model  of  the  female  breast, 
obtained  from  MRI,  within  which  we  have  assumed 
RI  variation  for  the  two  main  types,  namely,  adipose 
and  glandular  tissue.  Within  the  glandular  layer,  we 
have  inserted  an  abnormality  with  absorption  and 
scatter-only  contrast.  Using  the  modeled  NIR  data 
perturbed  by  noise,  we  have  reconstructed  images  of 
internal  absorption  and  reduced  scatter  and  have 
shown  the  effect  on  RI  variation  within  reconstructed 
images  for  the  cases  of  either  no  a  priori  or  exact 
information  on  the  RI  values. 

2.  Theory 

A.  Forward  Model 

It  is  generally  accepted  that  if  the  magnitude  of  the 
isotropic  fluence  rate  within  tissue  is  significantly 
larger  than  the  directional  flux  magnitude,  the  light 
field  is  diffuse,  which  occurs  when  the  scattering  in¬ 
teraction  dominates  over  absorption  in  a  region  of 
interest.  Mathematically,  this  assumption  allows  a 


simplification  of  the  Boltzmann  transport  equation, 
by  converting  the  description  of  an  anisotropic  light 
field  into  a  diffusion  equation  approximation.  The 
diffusion  approximation  in  the  frequency  domain  is 
given  by 


V  •  K(r)  V  J>(r,  go)  + 


^ a(r )  + 


i  oo 


cmir)_ 

X  0(r,  cd)  =  q0{r,  w),  (1) 


where  \ ±a  and  |jig'  are  absorption  and  reduced  scatter¬ 
ing,  respectively;  qQ{r,  go)  is  an  isotropic  source; 
J>(r,  go)  is  the  photon  fluence  rate  at  position  r; 


1 

K  “  3(p,a  +  p,s') 

is  the  diffusion  coefficient;  and  cm(r )  is  the  speed  of 
light  in  the  medium  at  any  point,  defined  by  c/n(r), 
where  n{r)  is  the  index  of  refraction  at  the  same  point 
and  c  is  the  speed  of  light  in  a  vacuum. 

The  best  description  of  the  air-tissue  boundary,  is 
reported  through  an  index-mismatched  type  III  con¬ 
dition,  in  which  the  fluence  at  the  edge  of  the  tissue 
exits  and  does  not  return.20  The  flux  leaving  the  ex¬ 
ternal  boundary  is  equal  to  the  fluence  rate  at  the 
boundary  weighted  by  a  factor  that  accounts  for  the 
internal  reflection  of  light  back  into  the  tissue.  This 
relation  is  described  in 


J>(^,  a))  +  2An  •  k(£)V  0(^,  w)  —  0,  (2) 

where  £  is  a  point  on  the  external  boundary  (dll!),  and 
A  depends  on  the  relative  RI  mismatch  between  tis¬ 
sue  ll!  and  air.  A  can  be  derived  from  Fresnel's  law, 


A  = 


2/(1  -  R0)  -  1  +  |  cos  0C|3 
1  -  |  cos  0C  | 3 


(3) 


where  0C  =  arcsin(ftAIR/ft!),  the  angle  at  which  total 
internal  reflection  occurs  for  photons  moving  from 
region  with  RI  n1  to  air  with  RI  nAIR,  and 


(ni/nAm  -  1) 
(nJriMB.  +  l)2' 


At  the  external  boundaries,  tiair  =  1,  the  RI  of  free 
space. 

At  interior  nodes,  which  lie  on  an  interface  between 
two  media  with  different  indices  of  refraction,  we 
apply  the  conditions  used  by  Schmitt  et  al.21  Taka- 
tani  and  Graham,22  and  Faris23: 


ii •D1V01(^,  (o)-n-D2V<b2(£,  (5) 


foi(£,  _  (ni\2 

<b2(£,  w)  \n2J  * 


(6) 
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Fig.  1.  (a)  Coronal  slice  of  an  MRI  of  a  female  subject  showing  adipose  and  glandular  tissue  layers  and  (b)  the  finite-element  model  mesh 

created  from  segmenting  the  MRI. 


These  equations  enforce  continuity  in  the  flux  across 
a  step  change  in  RI  (n),  and  establish  a  corresponding 
discontinuity  in  the  fluence  on  the  basis  of  the  two 
RIs  defining  the  regions  separating  the  boundary. 

B.  Inverse  Model 

The  data  is  represented  by  a  nonlinear  operator  y* 
=  F[|xa,  k],  where  y*  is  a  complex  vector  having  a 
real  and  imaginary  components,  which  are  mapped  to 
log  amplitude  and  phase  in  measurement.  Then  the 
image  reconstruction  method  seeks  a  solution  to 

(|la,  k)  =  arg  min^  J|  [y*  —  F([ia,  k)]||,  (7) 

where  ||.||  is  the  weighted  L2-norm,  representing  the 
square  root  of  the  sum  of  the  squared  elements.  The 
magnitude  of  this  is  sometimes  referred  to  as  the 
projection  error  and  provides  a  value  for  determining 
the  convergence  of  the  iterative  reconstruction  algo¬ 
rithm. 

The  finite-element  method  is  used  as  a  general  and 
flexible  method  for  solving  the  forward  problem  in 
arbitrary  geometries.11’24’25  In  the  inverse  problem,  in 
which  the  aim  is  to  recover  internal  optical  property 
distributions  from  boundary  measurements,  we  as¬ 
sume  that  fxa(r)  and  K(r)  are  expressed  in  a  basis  with 
a  limited  number  of  dimensions  (less  than  the  dimen¬ 
sion  of  the  finite-element  system  matrices).  A  number 
of  different  strategies  for  defining  reconstruction 
bases  are  possible;  in  this  paper  we  use  a  linear  pixel 
basis.  To  find  (|ia,  k)  in  Eq.  (7),  we  used  a  Levenberg- 
Marquardt  algorithm  in  which  is  repeatedly  solved 

a  =  JT(JJT+p/r1b,  (8) 

where  b  is  the  data  vector,  b  =  (y*  -  F[|xa,  k])t;  a  is 


the  solution  update  vector,  a  =  [8k;  8|xJ;  p  is  the 
regularization  factor;  and  J  is  the  Jacobian  (sensitiv¬ 
ity  or  weight)  matrix  for  the  DA  model  that  is  calcu¬ 
lated  with  the  Adjoint  method.26 

3.  Simulation  Methods 

An  MRI  of  a  female  subject,  Fig.  1(a),  was  used  to 
create  a  two-dimensional  mesh  Fig.  1(b).19  The  mesh 
contained  2306  nodes  corresponding  to  4042  linear 
triangular  elements  and  has  a  major  and  minor  axis 
diameter  of  101  and  82  mm,  respectively.  For  this 
subject,  experimental  data  was  collected  with  an 
MRI-NIR  system.19  We  combined  the  NIR  measure¬ 
ments  with  a  priori  anatomical  information  from  the 
MRI  (assuming  no  RI  variation)  to  estimate  the  op¬ 
tical  absorption  and  reduced  scatter  for  each  tissue 
type.  These  procedures  generated  estimates  of  the 
adipose  layer  as  having  an  absorption  of  0.003  mm-1 
and  a  reduced  scatter  of  0.95  mm-1,  whereas  the 
glandular  tissue  resulted  in  an  absorption  of 
0.006  mm-1  and  a  reduced  scatter  of  1.1  mm-1.  Al¬ 
ternatively,  assuming  a  single  tissue  type  within  the 
breast,  we  calculated  bulk  properties  for  the  breast  as 
having  an  absorption  of  0.0056  mm-1  and  a  reduced 
scatter  of  1.125  mm-1.  Assuming  an  RI  of  1.455  for 
the  adipose  tissue15  and  an  RI  of  1.2,  1.3,  or  1.4  for 
glandular  tissue,  we  calculated  data  with  the  pres¬ 
ence  of  an  anomaly  within  the  breast,  for  two  sepa¬ 
rate  cases: 

1 .  Assuming  a  homogenous  distribution  of  absorp¬ 
tion  and  reduced  scatter  of  0.0056  and  1.125  mm-1, 
respectively,  it  was  assumed  that  the  RI  varied  be¬ 
tween  each  layer  with  a  value  of  1.455  for  the  adipose 
tissue  and  1.2,  1.3,  or  1.4  for  glandular  tissue.  The 
modeled  anomaly,  as  shown  in  Fig.  2  (top  row)  had  a 
contrast  in  absorption  and  reduced  scatter  only  and 
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Fig.  2.  Images  [left  column,  absorption  (mm-1);  right  column, 
reduced  scatter  (mm-1)]  reconstructed  from  synthetic  data  gener¬ 
ated  with  a  model  with  RI  variation.  Top  row  shows  the  exact 
optical  properties;  the  adipose  and  glandular  zones  have  the  iden¬ 
tical  optical  properties  but  different  RI  of  1.455  and  1.4,  respec¬ 
tively.  Middle  row  shows  reconstructed  images  assuming  the 
correct  RI  (same  as  in  top  row)  distribution.  Bottom  row  shows 
reconstructed  images  assuming  a  homogenous  RI  of  1.455. 


was  assumed  to  have  the  same  RI  as  glandular  tis¬ 
sue. 

2.  The  adipose  layer  had  an  absorption  of 


0.003  mm-1  and  a  reduced  scatter  of  0.95  mm-1, 
whereas  the  glandular  tissue  had  an  absorption  of 
0.006  mm-1  and  a  reduced  scatter  of  1.1  mm-1.  The 
RI  varied  between  each  layer  with  a  value  of  1.455  for 
the  adipose  tissue  and  1.2,  1.3,  or  1.4  for  glandular 
tissue.  The  modeled  anomaly,  as  shown  in  Fig.  5 
below,  (top  row)  had  a  contrast  in  absorption  and 
reduced  scatter  only  and  was  assumed  to  have  the 
same  RI  as  glandular  tissue. 

Boundary  data  was  modeled  for  16  equally  spaced 
sources  and  detectors  arranged  on  the  external  pe¬ 
riphery  of  the  model,  with  added  Gaussian  noise  of 
1%  in  amplitude  and  1°  in  phase.  For  each  set  of  data, 
two  types  of  reconstruction  were  investigated:  (i)  cor¬ 
rect  (a  priori)  information  on  RI  variation  and  (ii)  a 
constant  RI  of  1.455.  In  all  cases  images  were  recon¬ 
structed  on  a  20  X  20  uniform  grid  basis  by  use  of  an 
initial  regularization  [Eq.  (8)]  of  p  =  10,  and  the 
number  of  iterations  was  continued  until  the  projec¬ 
tion  error  [Eq.  (7)]  changed  by  less  than  1%  from  the 
previous  iteration.1 

4.  Results 

In  each  of  the  cases  described  in  Section  3,  the  sim¬ 
ulated  data  were  calibrated  with  methods  described 
elsewhere127’28  to  provide  an  estimate  of  the  bulk 
optical  properties  of  absorption  and  reduced  scatter 
for  initiating  the  reconstruction  algorithm. 

Figure  2  shows  the  exact  (top  row)  and  the  recon¬ 
structed  images  for  the  model,  which  has  a  homoge¬ 
nous  absorption  and  reduced  scatter  of  0.0056  and 
1.125  mm-1,  respectively,  except  within  the  anomaly 
that  was  assigned  values  of  0.012  and  2  mm-1.  The  RI 
is  1.455  for  the  adipose  tissue  and  1.4  for  the  glan¬ 
dular  and  the  anomaly  tissues.  Images  were  recon¬ 
structed  assuming  either  correct  values  of  RI  for  each 
layer  (middle  row)  or  a  constant  RI  distribution  of 
1.455  (bottom  row).  In  each  case,  images  of  absorp¬ 
tion  and  reduced  scatter  were  recovered  simulta¬ 
neously  by  use  of  log  amplitude  and  phase  of  the  NIR 
data. 

Figures  3  and  4  present  exact  (top  row)  and  recon¬ 
structed  images  in  which  the  model  is  identical  to 
that  in  Fig.  2,  except  for  the  RI  of  the  glandular  and 
the  anomaly  tissues,  which  was  set  to  1.3  and  1.2, 
respectively. 

Figure  5  shows  exact  (top  row)  and  reconstructed 
images  in  which  the  adipose  layer  had  an  absorption 
of  0.003  mm-1  and  a  reduced  scatter  of  0.95  mm-1 
whereas  the  glandular  tissue  had  an  absorption  of 
0.006  mm-1  and  a  reduced  scatter  of  1.1  mm-1.  The 
anomaly  within  the  glandular  tissue  had  an  absorp¬ 
tion  of  0.012  mm-1  and  a  reduced  scatter  of  2  mm-1. 
The  RI  was  1.455  for  the  adipose  tissue  and  1.4  for  the 
glandular  and  the  anomaly  tissues.  As  in  Fig.  2,  im¬ 
ages  were  produced  assuming  either  the  correct  val¬ 
ues  for  the  RI  in  each  layer  (middle  row)  or  a 
homogeneous  RI  distribution  of  1.455  (bottom  row). 
Figures  6  and  7  show  analogous  results  to  those  in 
Figure  5,  except  for  the  RI  of  the  glandular  and  the 
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Fig.  3.  Same  as  Fig.  2  except  that  the  glandular  tissue  (and  the 
anomaly)  has  an  RI  of  1.3. 


Fig.  4.  Same  as  Fig.  2  except  that  the  glandular  tissue  (and  the 
anomaly)  has  an  RI  of  1.2. 


anomaly  tissues,  which  was  set  to  1.3  and  1.2,  respec¬ 
tively. 

5.  Discussion 

In  cases  in  which  the  tissue  is  optically  (absorption 
and  reduced  scatter)  homogenous  (except  for  an 
anomaly)  but  has  different  RIs  for  different  tissue 
types,  the  recovered  optical  properties  are  very  sim¬ 
ilar  to  the  expected  values,  yielding  a  peak  value  of 


the  anomaly  for  absorption  and  reduced  scatter  of 
approximately  0.0087  and  2.0  mm-1,  respectively, 
when  correct  a  priori  information  on  RI  distribution 
is  applied  throughout  the  model.  Assuming  that  the 
model  is  also  uniform  in  RI,  the  quality  of  the  recon¬ 
structed  images  varies,  depending  on  the  degree  of 
deviation  from  the  true  values  of  RI  as  expected.  For 
small  deviations  of  RI,  that  is,  when  the  glandular 
tissue  has  an  RI  of  1.455  but  is  assumed  to  be  1.4,  the 
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Fig.  5.  Images  [left  column,  absorption  (mm-1);  right  column, 
reduced  scatter  (mm-1)]  reconstructed  from  synthetic  data  gener¬ 
ated  with  a  model  with  RI  variation.  Top  row  shows  the  exact 
optical  properties;  the  adipose  and  glandular  tissue  have  different 
optical  properties  and  different  RI  of  1.455  and  1.4,  respectively. 
Middle  row  shows  reconstructed  images  assuming  the  correct  RI 
distribution.  Bottom  row  shows  reconstructed  images  assuming  a 
homogenous  RI  of  1.455. 


reconstructed  images  are  very  similar  to  the  case  in 
which  correct  RI  is  assumed,  Fig.  2.  As  the  deviation 
is  increased  (Figs.  3  and  4),  the  quality  of  the  recon¬ 


Fig.  6.  Same  as  Fig.  5  except  that  the  glandular  tissue  (and  the 
anomaly)  has  an  RI  of  1.3. 


structed  images  is  reduced.  For  example,  in  Fig.  4  (an 
assumed  value  of  1.2  for  glandular  tissue  RI),  the 
amount  of  background  noise  in  absorption  and  scat¬ 
ter  is  increased,  and  the  recovered  anomaly  is  blurred 
and  overly  smoothed.  Interestingly,  the  peak  value  of 
absorption  has  increased  to  0.01  mm-1. 

When  the  tissue  is  heterogeneous  in  both  optical 
(absorption  and  reduced  scatter)  properties  and  RI, 
correct  a  priori  information  on  the  RI  distribution 
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Fig.  7.  Same  as  Fig.  5  except  that  the  glandular  tissue  (and  the 
anomaly)  has  an  RI  of  1.2. 


throughout  the  model  generates  recovered  optical 
properties  that  are  very  similar  to  the  expected  val¬ 
ues,  yielding  a  peak  value  in  the  anomaly  for  absorp¬ 
tion  and  reduced  scatter  of  approximately  0.0086  and 
2.0  mm-1,  respectively.  If  the  reconstruction  assumes 
a  uniform  RI,  the  quality  of  the  reconstructed  images 
again  varies,  depending  on  the  level  of  deviation  from 
the  true  distribution  of  RI.  For  small  deviations,  the 


reconstructed  images  are  similar  to  the  case  in  which 
correct  RI  is  used  (Fig.  5).  As  the  deviation  of  the 
assumed  RI  increases  (Figs.  6  and  7),  the  quality  of 
the  reconstructed  images  is  reduced,  but  the  peak 
value  of  recovered  absorption  increases.  In  Fig.  7,  for 
example,  the  RI  of  the  glandular  tissue  is  1.2, 
whereas  the  assumed  value  for  the  homogenous  re¬ 
construction  is  1.455  (bottom  row).  As  in  the  homog¬ 
enous  case,  the  amount  of  background  noise  in 
absorption  and  scatter  is  increased;  however,  the  re¬ 
covered  anomaly  is  blurred  and  overly  smoothed,  and 
the  peak  value  of  absorption  has  increased  to 
0.01  mm1. 

Figure  8  shows  the  sensitivity  map  (Jacobian  ma¬ 
trix)  of  log  amplitude  and  absorption  for  a  single 
source  and  detector,  where  the  RI  distribution  is  as¬ 
sumed  to  be  either  homogeneous  at  1.455  or  hetero¬ 
geneous  with  an  adipose  layer  of  RI  =  1.455  and  a 
glandular  layer  of  RI  =  1.2.  From  this  figure,  it  evi¬ 
dent  that  a  decrease  in  RI  of  the  glandular  tissue 
reduces  the  sensitivity  deep  within  the  breast.  This 
reduction  in  sensitivity  may  explain  the  increased 
errors  observed  when  a  homogenous  assumption  of 
RI  is  used  for  image  reconstruction.  A  lower  sensitiv¬ 
ity  within  the  glandular  region  will  give  rise  to  overly 
smoothed  reconstructed  images.  However,  note  that 
an  assumption  for  the  glandular  tissue  RI  of  1.2  is 
probably  not  realistic  and  is  only  included  here  for 
bracketing  the  effects  on  image  reconstruction. 

The  recovery  of  the  reduced  scatter  images  seems 
to  have  been  less  affected  by  the  variation  of  the  RI. 
Although  only  the  sensitivity  maps  of  log  amplitude 
and  absorption  are  only  shown,  similar  trend  is  seen 
for  other  available  data  types  and  optical  properties 
(log  amplitude  or  phase  versus  absorption  and  scat¬ 
ter).  However,  it  is  important  to  note  that  the  mag¬ 
nitude  of  the  reduced  sensitivity  for  scatter  and  log 
amplitude  or  phase  (not  shown)  is  much  smaller  than 
that  seen  for  log  amplitude  and  absorption. 


6.  Conclusion 

The  effect  of  discrete  RI  changes  within  breast  tissue 
upon  NIR  image  reconstruction  has  been  investi¬ 
gated.  Previous  studies11  reported  our  results  in  mod¬ 
eling  the  effect  of  RI  on  the  forward  model.  In  our 
current  results,  the  initiative  was  advanced  to  inves¬ 
tigate  the  effect  of  RI  variation  on  NIR  image  recon¬ 
struction,  either  by  assumption  of  correct  knowledge 
of  the  RI  of  each  tissue  or  by  application  of  a  homog¬ 
enous  value  throughout  the  model. 

Synthetic  measurement  data  have  been  generated 
from  a  realistic  two-dimensional  breast  MRI,  exhib¬ 
iting  two  distinct  layers  of  adipose  and  glandular 
tissue.  Published  values  exist  for  the  RI  of  adipose 
tissue  but  not  for  the  glandular  tissue.  However,  it  is 
generally  accepted  that  glandular  tissue  has  a  lower 
RI,  and  we  have  assumed  a  value  of  1.4.16  Nonethe¬ 
less,  to  illustrate  the  effects  of  RI  variation  with  NIR 
image  reconstruction,  we  have  generated  data  as¬ 
suming  an  RI  of  glandular  tissue  of  1.4,  1.3,  or  1.2. 
Furthermore,  we  have  also  modeled  two  distinct 
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Fig.  8.  Sensitivity  map  (30  contour  lines)  of  log  amplitude  to  absorption  changes  for  a  single  source  and  detector  combination  for  (a)  a 
model  of  homogenous  RI  of  1.455  having  a  maximum  value  of  1.4  X  10-10  and  a  minimum  value  of  - 1.72  and  (b)  a  heterogeneous  RI  of 
adipose,  1.455,  and  glandular  tissue,  1.2,  having  a  maximum  value  of  1.3  X  10”10  and  a  minimum  value  of  -1.78. 


cases  in  which  each  layer  had  either  the  same  or 
different  absorption  and  reduced  scatter  properties, 
allowing  us  to  separate  the  effects  of  absorption  and 
reduced  scatter  from  RI  variation.  In  either  case,  we 
modeled  an  anomaly  deep  within  the  breast  and  re¬ 
constructed  images  assuming  either  correct  informa¬ 
tion  on  the  RI  or  simply  assuming  it  to  be 
homogenous. 

The  reconstructed  images  show  that  providing  the 
RI  of  the  glandular  tissue,  when  it  is  not  far  from  the 
value  of  adipose  tissue,  has  little  effect  on  the  quali¬ 
tative  and  quantitative  accuracy  of  the  results.  For 
the  layered  model  (Figs.  5-7)  the  difference  in  con¬ 
trast  of  the  reconstructed  anomaly  (between  use  of 
the  correct  RI  or  a  homogenous  value)  is  3.5%  when 
the  actual  glandular  tissue  RI  is  1.4  (but  assumed  to 
be  1.455)  and  increases  to  16%  when  the  actual  glan¬ 
dular  tissue  RI  is  1.2.  If  the  RI  variation  is  not  mod¬ 
eled,  the  background  noise  in  the  reconstructed 
images  increases,  and  the  reconstructed  anomaly  ex¬ 
hibits  a  more  blurred  character,  which  can  be  as  large 
as  200%  in  the  absorption  images.  This  increase  in 
noise  and  in  blurriness  can  be  understood  by  exami¬ 
nation  of  the  sensitivity  functions  in  Fig.  8,  which 
show  that  as  the  RI  of  the  glandular  tissue  decreases 
with  respect  to  adipose  tissue,  the  sensitivity  within 
the  glandular  tissue  also  decreases.  Another  interest¬ 
ing  and  very  important  point  to  note  here  is  that 
when  the  RI  distribution  is  not  correctly  modeled 
within  the  reconstruction  algorithm,  the  effects  of  the 
data  mismatched  owing  to  the  RI  are  exhibiting 
themselves  within  the  reconstructed  absorption  im¬ 
ages,  clearly  seen  in  Figs.  3  and  4.  Also,  images  were 
reconstructed  (not  shown)  assuming  a  homogeneous 
RI  of  1.3  and  1.2  in  all  of  the  presented  results  (de¬ 
pending  on  the  RI  of  the  glandular  tissue),  and  sim¬ 


ilar  results  were  found,  as  presented,  when  a 
homogeneous  RI  of  1.455  was  assumed. 

Our  research  has  considered  only  the  reconstruc¬ 
tion  of  absorption  and  reduced  scatter  with  data  gen¬ 
erated  by  a  model  with  RI  variation.  Assuming  the  RI 
of  the  whole  breast  is  similar,  for  example,  the  RI  is 
1.455  for  adipose  and  1.4  for  glandular  tissue,  recon¬ 
structed  images  of  absorption  and  scatter  can  be  ob¬ 
tained  that  ignore  the  effect  of  RI  with  modest 
degradation  in  the  recovery  of  information  about  an 
abnormality.  However,  it  is  also  important  to  note 
that  this  analysis  was  completed  under  the  assump¬ 
tion  that  the  abnormality  (i.e.,  tumor)  has  the  same 
RI  as  its  background  tissue  and  is  most  typically 
located  in  the  fibroglandular  tissue.  Although,  a 
study  by  Jiang  and  Xu12  has  reported  large  variations 
in  the  RI  of  tumors  relative  to  the  background,  little 
other  pathologic  or  in  vivo  data  exists  on  the  RI  of 
glandular  tissue  or  the  various  types  of  tumor.  Fur¬ 
ther  studies  are  needed  to  establish  the  RI  variation 
for  different  types  of  tumor  and  to  investigate  how 
such  variations  might  alter  NIR  tomography. 
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Spectrally  constrained  chromophore  and 
scattering  near-infrared  tomography  provides 
quantitative  and  robust  reconstruction 


Subhadra  Srinivasan,  Brian  W.  Pogue,  Shudong  Jiang, 
Hamid  Dehghani,  and  Keith  D.  Paulsen 


A  multispectral  direct  chromophore  and  scattering  reconstruction  technique  has  been  implemented  for 
near-infrared  frequency-domain  tomography  in  recovering  images  of  total  hemoglobin,  oxygen  satura¬ 
tion,  water,  and  scatter  parameters.  The  method  applies  the  spectral  constraint  of  the  chromophores  and 
scattering  spectra  directly  in  the  reconstruction  algorithm,  thereby  reducing  the  parameter  space  of  the 
inversion  process.  This  new  method  was  validated  by  use  of  simulated  and  experimental  data,  and  results 
show  better  robustness  and  stability  in  the  presence  of  higher  levels  of  noise.  The  method  suppresses 
artifacts,  especially  those  significant  in  water  and  scatter  power  images,  and  reduces  cross  talk  between 
chromophore  and  scatter  parameters.  Variation  in  scattering  was  followed  by  this  spectral  approach 
successfully  in  experimental  data  from  90-mm-diameter  cylindrical  phantoms,  and  results  show  linear 
variation  in  scatter  amplitude  and  reduced  scattering  coefficient  (|xs'),  with  total  hemoglobin,  oxygen 
saturation,  and  water  remaining  constant  and  quantitatively  accurate.  Similar  experiments  were  carried 
out  for  varying  oxygen  saturation  and  total  hemoglobin.  Accurate  quantification  was  obtained  with  a 
mean  error  of  7.7%  for  oxygen  saturation  and  6.2%  for  total  hemoglobin,  with  minimal  cross  talk  between 
different  parameters.  ©  2005  Optical  Society  of  America 
OCIS  codes:  170.3010,  170.6960. 


1.  Introduction 

Near-infrared  (NIR)  tomography  can  be  used  to  char¬ 
acterize  malignant  and  normal  tissue  based  on  the 
high-contrast  available  from  heme  in  the  blood,  lead¬ 
ing  to  images  that  are  related  to  intrinsic  pathophys¬ 
iologic  processes  such  as  angiogenesis  and  hypoxia. 
Absorption-based  parameters  can  be  recovered  such 
as  total  hemoglobin  in  the  tissue,  hemoglobin  oxygen 
saturation,  and  water  fraction.  It  is  also  possible  to 
estimate  elastic  scattering  images  that  may  provide 
information  about  the  composition  of  the  tissue.  In 
vivo  studies  have  demonstrated  levels  of  hemoglobin 
in  tumors  over  twice  that  in  normal  breast,1’2  and 
lower  levels  of  oxygen  saturation  have  been  found  in 
malignancies2’3;  however,  one  of  the  current  chal- 
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lenges  is  to  optimize  the  quantitative  accuracy  with 
which  these  parameters  can  be  determined.  The 
quantification  of  chromophores  and  scattering  pa¬ 
rameters  relies  upon  the  spectral  decomposition  of 
the  images  acquired  at  a  sparse  number  of  discrete 
wavelengths  instead  of  a  complete  spectrum.  This 
sparse  spectral  sampling  coupled  with  an  image  re¬ 
construction  process  that  is  ill  posed,  tends  to  amplify 
errors  in  quantifying  the  spatially  resolved  parame¬ 
ters  of  the  tissue.  In  this  study  a  spectrally  con¬ 
strained  approach  to  image  reconstruction  is 
introduced,  which  follows  the  recent  pioneering  de¬ 
velopments  proposed  by  Corlu  et  al.4  and  Li  et  al .,5 
who  showed  that  incorporation  of  spectral  informa¬ 
tion  into  the  reconstruction  process  improves  the 
uniqueness  of  the  image  formation  by  using 
continuous-wave  data.  In  this  paper  the  addition  of 
phase  information  and  the  improved  accuracy  in  fit¬ 
ting  water  and  scattering  power  are  specifically  ex¬ 
amined  by  use  of  this  spectrally  constrained 
approach.  The  improvement  in  quantification  of  wa¬ 
ter  and  scattering  has  dramatic  implications  in  terms 
of  the  value  of  these  particular  parameters  in  breast 
imaging.  The  potential  to  reduce  cross  talk  between 
chromophores  is  also  important. 
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In  earlier  research  the  absorption  and  scattering 
coefficients  were  recovered  from  boundary  measure¬ 
ments  of  amplitude  or  phase  or  both  on  the  object  of 
interest,  such  as  the  breast  or  brain,  by  means  of 
computational  models  in  which  both  analytical  and 
numerical  approaches  were  used  by  different  re¬ 
search  groups.6-8  After  recovery  of  these  optical  co¬ 
efficients,  a  spectral  fitting  to  known  absorption 
signatures  of  oxyhemoglobin,  deoxyhemoglobin,  and 
water  is  used  to  obtain  the  concentrations  of  these 
absorbing  chromophores.  Similarly,  the  reduced  scat¬ 
tering  coefficients  at  separate  wavelengths  were  fit  to 
yield  the  scatter  amplitude  (a)  and  scatter  power  (6), 
which  are  related  to  the  structure  of  the  tissue  in 
terms  of  scatterer  size  and  density.  In  this  paper  a 
modified  reconstruction  approach  is  used,  which  im¬ 
plements  the  possible  spectral  shapes  of  the  chro- 
mophore  and  scattering  models  into  the  image 
formation  process,  thereby  adding  a  spectral  con¬ 
straint  into  the  reconstruction.  The  chromophore  con¬ 
centrations  and  scatter  parameters  are  estimated 
directly  by  incorporating  the  known  Beer’s  law  atten¬ 
uation  relation  and  Mie  scattering  behavior  as  con¬ 
straints.  This  type  of  parameter  reduction  has  been 
applied  in  electrical  impedance  tomography  in  which 
Brandstatter  et  al.9  showed  that  by  using  multifre¬ 
quency  data  and  by  assuming  a  frequency  depen¬ 
dence,  one  can  reduce  the  ill-posed  nature  of  the 
problem  and  make  the  reconstruction  more  robust  to 
noise  in  data.  A  similar  application  in  microwave 
image  reconstruction10  provides  evidence  of  reduced 
artifacts  in  the  images  as  a  result  of  coupling  mea¬ 
surements  from  different  frequencies.  Corlu  et  al.4 
implemented  this  approach  by  using  continuous- 
wave  (cw)  measurements  to  find  the  optimal  four 
wavelengths  that  reduce  the  cross  talk  between  ab¬ 
sorption  and  scatter  parameters.  Their  results  from 
simulations  are  encouraging  and  are  based  on  the 
assumption  that  all  change  in  scattering  is  due  to  the 
scatter  amplitude  with  the  scatter  power  kept  con¬ 
stant.  A  similar  approach  to  cw  data  has  been  imple¬ 
mented  by  Li  et  al.5  used  two  of  three  wavelengths 
under  the  assumption  that  there  is  no  scattering  per¬ 
turbation.  They  have  applied  this  method  to  find 
chromophore  concentrations  directly  and  have  shown 
in  simulated  and  experimental  data  that  their  tech¬ 
nique  results  in  reduced  image  artifacts  and  param¬ 
eter  cross  talk. 

In  the  current  study  this  overall  approach  is  ex¬ 
tended  to  the  application  of  frequency-domain  data, 
using  six  wavelengths.  The  method  is  evaluated  with 
experimental  data  following  individual  variation  of 
oxygen  saturation,  hemoglobin,  and  scattering  pa¬ 
rameters.  A  finite-element  model  of  the  diffusion 
equation  is  used,  and  the  algorithm  reconstructs  im¬ 
ages  for  five  parameters:  oxyhemoglobin,  deoxyhemo¬ 
globin,  water  fraction,  scatter  amplitude,  and  scatter 
power,  with  no  assumptions  on  the  scatter  amplitude 
or  power.  The  results  show  that  the  new  technique  is 
more  robust  to  noise  in  measurements  than  the  con¬ 
ventional  method.  In  addition,  the  spectral  constraint 
reduces  the  noise  in  the  recovered  chromophore  con¬ 


centrations,  especially  in  the  water  and  scattering 
images,  and  the  reconstructions  from  the  experimen¬ 
tal  data  show  quantitatively  accurate  results. 

2.  Materials  and  Methods 

A.  Instrumentation 

The  NIR  frequency-domain  system  for  breast  imag¬ 
ing  has  been  described  in  detail  in  previous  papers.11 
Briefly,  it  consists  of  optical  fibers  placed  in  three 
planes  in  a  circular  geometry.  Each  plane  has  16 
source- detector  positions,  and  intensity  light  modu¬ 
lated  at  100  MHz  is  used  at  six  different  wavelengths 
in  the  range  660-850  nm.  The  signals  are  detected  by 
high-gain  photomultiplier  tubes,  and  the  electrical 
signals  are  passed  through  rf  mixer  circuits  to  het¬ 
erodyne  down  to  a  500-Hz  offset  frequency.  The  sig¬ 
nal  amplitude  and  phase  are  calibrated  to 
compensate  for  system  offsets  by  matching  measured 
data  from  homogeneous  phantoms12  to  simulated  re¬ 
sults  from  the  finite-element  model.  When  optimized, 
the  calibrated  data  have  less  than  1%  offset  from 
simulated  values  and  provide  a  highly  stable  data  set 
from  which  to  reconstruct  absorption  and  scattering 
coefficient  images. 

B.  Reconstruction  without  Spectral  Constraints 

Under  the  assumption  that  breast  tissue  is  a  scatter- 
dominated  medium,  the  diffusion  approximation  to 
the  radiative  transfer  equation8’13  was  used  to  model 
the  propagation  of  light  at  large  distances  from  the 
source  location.  This  is  given  by 

/  m\ 

-  V  •  K (r)  V  0(r,  go)  +  ( fxa(r)  +  —  )<l>(r,  go)  =  q0(r,  go), 

(1) 

where  0(r,  go)  is  the  isotropic  fluence  at  modulation 
frequency  go  and  position  r,  K(r)  is  the  diffusion  coef¬ 
ficient,  | xa(r)  is  the  absorption  coefficient,  c  is  the 
speed  of  light  in  the  medium,  and  q0(r ,  go)  is  an  iso¬ 
tropic  source.  The  diffusion  coefficient  can  be  written 
as 

1 

K(r)  =  3[jjLa(r)+,V(r)]’  (2) 

where  |xs'  is  the  reduced  scattering  coefficient. 

A  finite-element-based-model  solution  to  this  equa¬ 
tion  was  developed  and  was  described  and  validated 
in  Refs.  14  and  15.  Briefly,  the  forward  problem  in¬ 
volves  solving  Eq.  (1)  for  fluence,  given  an  initial 
distribution  of  the  optical  properties,  with  the  appro¬ 
priate  boundary  conditions  applied.  The  image  recon¬ 
struction  uses  a  Newton-Raphson  minimization  that 
iteratively  updates  the  optical  property  parameters 
based  on  a  least-squares-error  norm  given  by 
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where  c|);meas  is  the  measured  data  and  c| );cal  is  the  data 
calculated  for  an  initial  distribution  of  the  properties, 
using  the  forward  solver.  Here  the  measurements  are 
assumed  to  be  independent  of  each  other.  Using  a 
Taylor’s  series  approximation  for  the  solution  at  some 
close  distance  from  the  boundary  data  (4>cal)  for  the 
initial  distribution  and  ignoring  the  higher-order 
terms,  the  update  in  the  optical  properties  is  related 
to  the  difference  between  the  measured  and  the  cal¬ 
culated  data  as 

d<|>  =  $d|x,  (4) 

where  d§  refers  to  the  change  in  boundary  data;  $  is 
the  Jacobian,  the  matrix  containing  the  sensitivity  of 
the  boundary  data  to  a  change  in  optical  properties, 
$  —  C\Va;$K];  and  dp,  is  the  update  in  the  optical 
properties  given  by  dp,  =  [dp,a;di<]. 

The  reconstruction  is  sensitive  to  the  initial  esti¬ 
mate  of  the  parameters,  which  are  given  by  a  homo¬ 
geneous  prefitting  algorithm  based  on  the  analytical 
solution  for  infinite  medium.12  The  matrix  $,  being  ill 
conditioned,  requires  that  the  inverse  problem  in  Eq. 
(4)  be  solved  with  the  application  of  a  Levenberg- 
Marquardt  regularization  scheme1617  for  stabiliza¬ 
tion.  The  stopping  criterion  for  this  reconstruction 
was  chosen  to  be  when  the  y2  error  in  Eq.  (3),  known 
as  the  projection  error,  changes  by  less  than  2%  be¬ 
tween  successive  iterations. 

Previously18  the  optical  properties  at  each  wave¬ 
length  were  obtained,  and  then  the  calculation  of  the 
chromophore  concentrations  was  performed  with  a 
constrained  least-squares  fit  to  the  Beer’s  law  rela¬ 
tion 

|x0  =  [e]c,  (5) 


|xs'=a\-6. 

Equation  (6)  was  used  to  derive  the  scatter  amplitude 
(a)  and  the  scatter  power  (6),  with  wavelength  in 
micrometers.  The  coefficient  p,s'  has  units  of  inverse 
millimeters,  and  b  is  dimensionless  so  that  a  has 
units  given  by  l(T36(mm)6-1.  Both  the  scattering 
power  and  the  amplitude  depend  on  the  scattering 
center  size  and  number  density  and  may  reflect  vari¬ 
ations  in  breast  structural  composition  due  to  differ¬ 
ent  cellular,  organelle,  and  structural  sizes  and 
densities  for  fatty  and  glandular  tissue.  Typically, 
large  scatterers  have  lower  b  and  a  values,  whereas 
small  scatterers  have  higher  b  and  a  coefficients.19’20 
Although  scatterers  in  tissue  are  not  necessarily  ho¬ 
mogeneous  spheres,  as  assumed  in  Eq.  (6),  studies  on 
red  blood  cells  and  yeast  have  shown  this  to  be  a 
reasonable  approximation  since  measurements  of  ps' 
in  these  cells  show  trends  similar  to  Mie  theory.21’22 

C.  Spectrally  Constrained  Chromophore  and  Scattering 
Reconstruction 

Instead  of  estimating  the  optical  properties  at  each 
wavelength  and  then  spectrally  deconvolving  the 
chromophore  concentrations  [Eqs.  (5)  and  (6)],  one 
can  incorporate  these  constraints  into  the  reconstruc¬ 
tion  to  directly  determine  oxyhemoglobin,  deoxyhe¬ 
moglobin,  water,  scatter  amplitude,  and  scatter 
power,  thus  reducing  the  parameter  space  from  12 
images  (p,a  and  p,s'  at  6  wavelengths)  to  5  parametric 
images.  Assuming  that  we  know  £5^  =  d4>/dp,  and 
$K  =  d4>/dK,  as  calculated  by  the  previous  method 
(reconstruction  without  spectral  priors),  in  the  new 
approach  the  measurements  at  all  wavelengths  are 
coupled  together,  and  the  relations  in  Eqs.  (5)  and  (6) 
are  combined  to  create  a  new  set  of  relations,  which 
for  each  wavelength  is  represented  by 


where  e  is  the  molar  absorption  spectra  of  the  absorb¬ 
ing  chromophores  and  c  is  the  concentration  of  these 
chromophores.  Oxyhemoglobin  (Hb02),  deoxyhemo¬ 
globin  (Hb),  and  water  are  assumed  to  be  the  main 
absorbers,  and  their  molar  absorption  spectra  were 
obtained  experimentally  in  our  instrument.  This  ap¬ 
proach  of  using  values  estimated  by  the  system  com¬ 
pensates  for  any  offsets  from  the  theoretical  values, 
yet  there  was  little  difference  between  our  experi¬ 
mental  and  theoretical  estimates  of  molar  absorption 
coefficients.  By  fitting  for  the  concentrations,  we  cal¬ 
culate  total  hemoglobin  as  HbT  =  Hb02  +  Hb  [in 
micromolar  (p,M)]  and  oxygen  saturation  as  S02 
=  Hb02/HbT  X  100  (in  percent);  the  contribution 
from  other  chromophores  such  as  lipids  has  been 
found  to  be  negligible,  because  the  wavelengths  used 
here  were  limited  to  less  than  850  nm  where  lipid  is 
a  weak  absorber.  The  constraints  on  the  fitting  pro¬ 
cess  were  for  HbT  to  be  below  100  p,M,  oxygen  satu¬ 
ration  to  be  nonnegative  and  with  an  upper  bound  of 
100%,  and  water  to  be  in  the  range  [0%,  100%]. 

Similarly,  the  p,s'  spectrum  of  tissue  has  been 
shown  to  fit  well  to  an  empirical  approximation  to 
Mie  scattering  theory,19  20  given  by 


=  &,x  d  c  +  $0>x  d  a  +  !3&,x  d  b, 


(7) 


where 


dc 


<94>  dp, 
d[X  dc  ’ 


for  each  chromophore  (c)  in  the  model.  From  Eq.  (5) 
we  get  dp,  =  sdc,  so  that,  substituting  for  dp,/dc, 


cv  = 

Oc,X 


d  (f) 
dc 


\ 


d  (J) 

<8>  (Bxc1,  c2,  c3)  =  x  ®  (Sxa’ c2,  c3), 


(8) 


where  ®  refers  to  the  Kronecker  tensor  product. 
Similarly, 


r\Sa,\  = 


d  (J) 
da 


d(|)  d  K 

dK  da 


(9) 


I860 


APPLIED  OPTICS  /  Vol.  44,  No.  10  /  1  April  2005 


Rewriting 


3k 

(dK\ 

3a 

l  da  ) 

and  knowing  that 


1 

3(|xa  +  |V)’ 


we  get 


dK 

3|V 

da 


1 

3 

=  X” 


-1 

(g,a  +  iO2 


=  3  (-9k2) 


=  -3k2 


Substituting  these  expressions  in  Eq.  (9)  leads  to 
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Similarly,  for  the  scatter  power 
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Knowing  that  3  In  p,s'  =  (l/|x/)  d|x/^d|V/d  hi  |Ag' 
=  |jls'  and  from  Eq.  (6),  In  p,s'  =  In  a  -  b  In  A,  then  it 
is  found  that  3  In  p^'/dfe  =  -  In  X.  Substituting  these 
relations  produces 

$&,x  =  $k(  — 3k2)(p,s')(— In  \)  |  x.  (12) 

The  overall  system  of  equations  is  assembled  by 
substituting  the  relations  from  Eqs.  (8),  (10),  and  (12) 
into  Eq.  (7): 
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The  size  of  the  left-hand  vector  is  equal  to  the  number 
of  wavelengths  multiplied  by  the  number  of  measure¬ 
ments  per  wavelength,  and  the  size  of  the  right-hand 
side  vector  is  equal  to  the  number  of  chromophores 
and  scatter  parameters  multiplied  by  the  number  of 
nodes  for  updating  each  parameter  in  the  mesh.  The 
individual  Jacobian  matrices  for  each  chromophore 
were  computed  with  a  dual-mesh  technique,23  on  a 
2000-node  mesh  for  forward  diffusion  calculations 
and  a  425-node  mesh  was  used  for  image  reconstruc¬ 


tion.  Equation  (13)  was  implemented  by  building  the 
new  Jacobian  (5  X  425  by  480  X  6),  and  the  same 
Levenberg-Marquardt  regularization  scheme  was 
applied.  The  computation  time  was  approximately  25 
min  for  typically  5-7  iterations,  with  the  measure  of 
convergence  being  when  the  projection  error  was  less 
than  2%  of  the  previous  iteration  value.  Additional 
constraints  based  on  the  physiologically  possible  val¬ 
ues  for  the  parameters  were  applied  at  each  iteration 
so  that  HbT  does  not  exceed  100  p,M  (based  on  typical 
concentrations  found  in  the  breast),  that  oxygen  sat¬ 
uration  is  in  the  range  [0%,  100%] ,  and  that  water  is 
in  the  range  [0%,  100%].  The  scatter  amplitude  is 
bounded  in  the  range  [0.5,  2.0]  in  units  of 
10~36(mm)6_1,  and  the  scatter  power  is  in  the  range 
[0,  2]  based  on  previous  studies,24  so  that  together 
they  cover  the  possible  range  for  the  reduced  scatter¬ 
ing  coefficient.  This  range  is  [0.5,3.25]  mnT1  for 
785  nm.  The  approach  can  easily  be  extended  to  ad¬ 
ditional  wavelengths  without  any  computational  ex¬ 
pense  in  the  inversion  process  since  the  size  of  the 
new  Hessian  from  Eq.  (13)  depends  on  the  number  of 
nodes  and  not  on  the  number  of  measurements,  al¬ 
though  the  number  of  wavelengths  will  influence  the 
calculation  of  individual  Jacobian  matrices.  The  al¬ 
gorithm  typically  converges  after  a  lower  number  of 
iterations  than  the  conventional  method  does,  and  no 
spatial  filtering  was  necessary  since  the  noise  in  the 
images  is  already  damped  by  the  spectral  constraints. 

3.  Results 

Subsections  3.A  and  3.B  show  the  improvement  of  the 
spectral  technique  over  the  conventional  method  by 
quantitative  assessment  of  the  mean  and  standard 
deviation  from  recovered  images  by  use  of  simulated 
and  experimental  data.  The  results  shown  in  Subsec¬ 
tions  3.C,  3.D,  and  3.E  focus  on  validating  the  accu¬ 
racy  of  the  spectral  technique  in  following  the 
variation  of  scattering,  oxygen  saturation  (Hill 
curve),  and  total  hemoglobin  individually  by  use  of 
the  appropriate  experiments. 

A.  Effect  of  Noise  in  Amplitude  and  Phase 
Measurements 

It  is  expected  that  the  spectrally  constrained  nature 
of  this  new  algorithm  will  make  the  reconstruction 
more  immune  to  noise  in  measurements  as  compared 
with  the  conventional  method,  since  all  data  are  cou¬ 
pled.  More  noise  in  the  data  is  typically  observed  at 
661  nm  for  tissues  containing  lower  oxygenation  (due 
to  high  absorption  of  deoxyhemoglobin)  and  849  nm 
(due  to  water  absorption)  and  in  higher-scattering 
cases.  This  behavior  is  taken  into  consideration  in  the 
algorithm  by  the  a  priori  spectral  information.  To  test 
the  hypothesis  of  reduced  noise  sensitivity  of  the 
spectral  method,  we  simulated  amplitude  and  phase 
data  at  six  wavelengths  by  using  the  finite-element 
model  for  a  homogeneous  phantom  of  diameter 
86  mm  with  concentrations  of  30  jxM  Hfc>02, 
30-jxM  Hb,  60%  water,  scatter  amplitude  of  1  [units 
of  10-36(mm)6-1],  and  scatter  power  of  1.  This  yields 
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Fig.  1.  Recovered  mean  values  with  standard  deviation  error  bars  are  shown  for  (a)  hemoglobin,  (b)  oxygen  saturation  (S02),  (c)  water 
(in  percent),  (d)  scatter  amplitude,  and  (e)  scatter  power.  These  were  estimated  from  the  interior  of  a  homogeneous  field  reconstructed  with 
different  levels  of  noise  in  the  original  data.  Values  for  the  new  spectrally  constrained  reconstruction  are  shown  alongside  results  from  the 
conventional  approach. 


total  hemoglobin  of  60  p,M,  S02  =  50%,  and  jxs'  at 
785  nm  =  1.27  mm-1,  which  are  concentrations  typ¬ 
ically  found  in  the  breast.  Random  Gaussian  distrib¬ 
uted  noise  was  added  to  the  amplitude  and  phase 
measurements  in  a  systematic  manner  from  0.5%  to 
5%,  and  the  spectrally  constrained  reconstruction 
was  carried  out  on  the  data.  The  conventional  tech¬ 
nique  of  reconstructing  each  wavelength  separately 
and  then  carrying  out  the  spectral  fit  was  also  applied 
to  these  data  for  comparison.  The  mean  and  standard 
deviation  for  the  reconstructed  images  from  both 
techniques  for  each  of  the  parameters  are  plotted  in 


Fig.  1;  the  results  are  shown  for  the  cases  with  no 
noise,  1%  noise,  and  5%  noise. 

For  the  noiseless  data  reconstruction,  both  tech¬ 
niques  show  an  accurate  recovery  of  all  five  parame¬ 
ters  (mean  is  within  3%  of  the  true  value),  with  an 
average  standard  deviation  of  0.5%  of  the  mean  for 
the  spectral  method  and  3.7%  for  the  conventional 
technique.  For  the  1%  and  5%  noise  cases,  the  stan¬ 
dard  deviation  increases  as  expected;  however,  this 
increase  is  much  more  evident  in  the  images  from  the 
conventional  technique  than  from  the  spectral 
method.  The  results  from  the  spectral  method  do  not 
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differ  significantly  from  the  conventional  technique 
for  total  hemoglobin  and  oxygen  saturation  images; 
however,  the  noise  is  significantly  suppressed  in  the 
water  and  scatter  images  from  the  spectral  method. 
For  1%  noise  in  data,  the  mean  is  still  within  5%  of 
the  true  value  for  both  methods,  but  the  standard 
deviation  in  water  has  reduced  from  40%  in  the  con¬ 
ventional  method  to  12%  of  the  mean  for  the  spectral 
technique.  Even  in  the  5%  noise  case,  the  spectral 
method  shows  a  reasonable  recovery  of  mean  values 
for  the  parameters  (accurate  to  within  10%  on  an 
average),  with  a  15.3%  standard  deviation  (as  per¬ 
cent  of  the  mean).  This  shows  a  reduced  sensitivity  of 
the  reconstruction  to  higher  levels  of  noise  compared 
with  the  conventional  method. 

B.  Reduced  Standard  Deviation  in  Homogeneous 
Experimental  Data 

To  assess  the  mean  and  standard  deviation  from  ho¬ 
mogeneous  images  in  experimental  data,  we  collected 
measurements  on  a  liquid  tissue-simulating  phan¬ 
tom  within  a  plastic  circular  container  of  90-mm  di¬ 
ameter,  consisting  of  9.3-jxM  pig  blood  in  buffered 
saline  and  1%  Intralipid  concentration.  The  blood 
hematocrit  was  measured  before  the  experiment  with 
a  clinical  co-oximeter  that  showed  1%  of  the  pig  blood 
contained  9.3  jxM  of  hemoglobin  for  this  sample.  The 
expected  values  for  the  scatter  amplitude  and  scatter 
power  were  derived  from  the  work  of  van  Staveren  et 
al.19  Using  the  expression  given  by  van  Staveren  et 
al.  with  the  units  suitably  modified  produces  the  fol¬ 
lowing  equation:  p,s'  =  0.928X-14  -  0.16X-24.  Incor¬ 
porating  the  higher-order  term  into  the  scatter 
amplitude  factor  (since  the  amplitude  factor  of  the 
second  term  is  much  lower  than  that  of  the  first  term) 
by  assuming  that  p,s'  =  1  mm-1  at  800  nm,  the  scat¬ 
ter  amplitude  (a)  =  0.73  and  the  scatter  power  (b) 
=  1.4.  Water  and  oxygen  saturation  are  expected  to 
be  close  to  100%,  since  the  phantom  was  an  oxygen¬ 
ated  liquid  solution.  Both  the  spectral  and  conven¬ 
tional  techniques  were  applied  to  this  data,  and  the 
mean  and  standard  deviation  from  the  reconstructed 
images  are  plotted  in  Fig.  2(a),  along  with  the  ex¬ 
pected  values.  As  described  in  Section  2,  both  recon¬ 
structions  were  terminated  by  use  of  the  projection 
error  criterion,  and  reconstructed  parameters  15  mm 
from  the  edge  have  been  removed  from  the  calcula¬ 
tion  of  mean  and  standard  deviation  to  diminish  con¬ 
tribution  from  any  boundary  artifacts. 

Figure  2(a)  shows  the  reduced  standard  deviation  in 
the  images  obtained  from  the  spectral  method,  com¬ 
pared  with  the  conventional  technique.  The  mean  val¬ 
ues  for  the  parameters  are  accurate  to  within  6%,  on 
average,  for  the  spectral  scheme  and  to  within  11%  for 
the  conventional  method.  As  with  the  simulations,  no 
spatial  filtering  is  applied  to  the  spectral  reconstruc¬ 
tion,  whereas  the  conventional  method  uses  a  mean 
filter.  The  stopping  criterion  for  the  spectral  technique 
is  a  projection  error  change  of  less  than  2%  between 
iterations,  and  it  converges  in  seven  iterations.  For  the 
conventional  method,  the  equivalent  7th  iteration  at 


every  wavelength  was  used  to  obtain  images  based  on 
earlier  studies25  that  indicated  5-9  iterations  are  most 
suitable  for  experimental  data.  Both  methods  use  the 
same  initial  regularization  parameter  (equal  to  10  in 
this  study15).  The  main  improvement  here  was  the 
suppression  of  noise  in  the  water  and  scattering  im¬ 
ages  by  use  of  the  spectral  technique. 

Figure  2(b)  shows  the  reconstructed  images  from 
both  methods  along  with  a  cross  section  of  the  middle 
plane.  The  spatial  variation  in  the  cross  section  is  less 
in  the  spectral  technique,  and  some  of  the  boundary 
artifacts  in  the  hemoglobin,  water,  and  scatter  power 
images  are  reduced.  The  scatter  amplitude  image 
shows  some  artifacts  close  to  the  boundary,  which 
may  indicate  the  need  for  higher  regularization  for 
the  parameter.  The  hemoglobin  and  water  images 
from  the  conventional  technique  show  some  cross  talk 
between  the  images.  A  central  artifact  can  be  seen  in 
the  HbT  image  in  which  there  is  a  decrease  in  its 
value,  with  saturation  in  the  water  image  at  this 
same  region  (=  100%).  The  scatter  parameter  images, 
especially  scatter  power,  show  considerable  noise 
that  is  possibly  due  to  cross  talk  between  the  two 
scatter  parameters  and  between  deoxyhemoglobin 
and  scatter. 

C.  Scattering  Parameter  Validation 

Having  shown  that  spectral  reconstruction  is  supe¬ 
rior  to  the  conventional  method  of  reconstructing 
each  wavelength  separately  and  then  applying  the 
spectral  information  in  terms  of  reduced  sensitivity 
to  noise  in  the  data  and  suppression  of  artifacts  in  the 
images,  we  now  focus  the  following  sections  on  vali¬ 
dating  the  accuracy  of  the  spectral  reconstruction. 
One  of  the  key  advantages  of  the  spectral  method 
over  the  conventional  technique  is  the  reduction  of 
noise  in  the  water  and  scatter  parameters.  The  im¬ 
plementation  of  this  technique  on  frequency-domain 
measurements  allows  the  separation  of  absorption 
and  scatter,  and  this  along  with  multiwavelength  ca¬ 
pability  allows  a  modest  separation  of  the  scatter 
amplitude  and  scatter  power.  To  test  the  cross  talk 
and  noise  in  the  scatter  parameters  with  experimen¬ 
tal  measurements  and  to  follow  their  variation,  the 
Intralipid  concentration  was  varied  in  the  liquid 
phantom  solution  containing  1%  blood,  from  0.75%  to 
1.5%  in  steps  of  0.25%  (the  data  from  1%  Intralipid 
was  also  used  for  Fig.  2).  The  amplitude  and  phase 
measurements  were  taken  for  each  concentration, 
and  the  spectrally  constrained  reconstruction  was  ap¬ 
plied  to  the  data.  The  total  hemoglobin  was  constant 
through  the  varying  concentrations  of  Intralipid,  and 
the  saturation  for  both  water  and  oxygen  was  100% 
for  all  data  sets.  The  mean  value  along  with  standard 
deviation  from  the  images  are  plotted  for  scatter  am¬ 
plitude  and  scatter  power  in  Fig.  3(a),  and  Fig.  3(b) 
shows  the  average  p,s'  at  661  and  785  nm.  Figure  3(c) 
shows  the  total  hemoglobin  and  Fig.  3(d)  shows  the 
oxygen  saturation  and  water  content.  The  scatter  am¬ 
plitude  varies  linearly  with  concentration  and  shows 
more  variation  (range  0.6-1.25)  than  the  scatter 
power  (range  1.3-1.53).  Scatter  power  values  are 
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Fig.  2.  Recovered  mean  values  with  standard  deviation  error  bars  are  shown  in  (a)  from  reconstructed  data  with  a  90-mm-diameter  liquid 
phantom  containing  1%  Intralipid  with  9.3-|jlM  total  hemoglobin.  Values  for  the  spectrally  constrained  reconstruction  are  shown  alongside 
those  obtained  with  the  conventional  reconstruction  approach  and  the  true  theoretical  values,  (b)  Images  from  this  phantom  are  shown 
for  comparison  with  the  spectrally  constrained  reconstruction  (top  row),  the  conventional  reconstruction  (middle  row),  and  the  profile  plots 
from  the  midplane  of  these  images  (bottom  row). 


comparable  with  the  expected  value  of  1.4  from  van 
Staveren  et  al .,19  showing  a  mean  of  1.4  ±  0.1 
through  the  change  in  concentrations.  The  mean  jxs' 
at  785  nm  varies  linearly  (slope  of  ~1)  with  the 
change  in  percent  Intralipid,  and  the  value  approxi¬ 
mately  doubles  (0.89  versus  1.7  mm-1)  when  the  con¬ 
centration  doubles  from  0.75%  to  1.5%  Intralipid, 
which  is  encouraging.  The  reduced  scattering  coeffi¬ 
cient  at  661  nm  also  shows  a  similar  trend.  The  total 
hemoglobin  stays  constant  with  change  in  scattering 


with  a  mean  of  8.2  jxM  ±  0.8,  and  the  oxygen  satu¬ 
ration  shows  a  mean  value  of  99.3%  ±  1.2%,  close  to 
expected  value  of  100%.  Water  shows  an  average  of 
92.4%  ±  4.2%,  and  some  cross  talk  can  be  seen  be¬ 
tween  hemoglobin  and  water  at  the  higher  Intralipid 
concentrations  owing  to  the  high-scattering  medium. 

D.  Oxygen  Saturation  Validation 

Malignant  tumors  typically  have  lower  partial  pres¬ 
sure  values  for  oxygen  (P02)  owing  to  hypoxia,26  and 
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Fig.  3.  Recovered  mean  values  are  shown  from  a  series  of  phan¬ 
toms  in  which  the  scattering  coefficient  was  systematically  var¬ 
ied  through  concentrations  (Cone)  of  Intralipid  ranging  be¬ 
tween  0.75%  and  1.5%.  The  estimated  scattering  power  and 
amplitudes  are  shown  in  (a),  and  the  reduced  scattering  coeffi¬ 
cients  at  661  and  785  nm  are  shown  in  (b).  The  total  hemoglobin, 
which  did  not  vary,  is  shown  in  (c),  along  with  a  line  correspond¬ 
ing  to  the  theoretical  value.  In  (d)  the  oxygen  saturation  and 
water  values  are  shown,  which  also  did  not  vary.  Both  have 
theoretical  estimates  of  100%.  Error  bars  represent  the  standard 
deviation  of  all  pixels  within  the  interior  60  mm  of  the  region 
imaged. 


it  is  useful  to  verify  that  the  spectrally  constrained 
reconstruction  can  follow  these  hypoxic  conditions. 
The  relation  between  oxygen  saturation  and  the  vari¬ 
ation  in  partial  pressure  of  oxygen  has  a  well- 
characterized  behavior  given  by  the  Hill  curve.  This 
curve  was  obtained  previously  for  the  conventional 
method  and  published  in  Ref.  27.  Data  were  acquired 
by  use  of  a  phantom  solution  containing  1%  whole 
blood  and  1%  Intralipid  in  saline,  in  a  thin-walled 
plastic  container,  70-mm  in  diameter.  The  whole 
blood  (1%)  was  found  to  have  18-jxM  hemoglobin,  and 
the  oxygenation  of  the  solution  was  reduced  by  vary¬ 
ing  the  P02  values  from  150  to  0  mm  Hg  by  the  ad¬ 
dition  of  yeast.  The  P02  was  independently  measured 
by  means  of  a  chemical  microelectrode  after  calibra¬ 
tion  of  the  electrode  overnight  in  a  saline  solution.  By 
varying  the  P02  gradually,  using  a  small  amount  of 
yeast,  and  making  measurements  over  this  period  of 
time,  we  eventually  reduced  the  P02  to  zero  and  ob¬ 
tained  a  complete  set  of  data  over  the  required  range. 
The  spectral  reconstruction  of  this  data  gave  HbT, 
S02,  water,  scatter  amplitude,  and  scatter  power  im¬ 
ages,  from  which  the  mean  and  standard  deviation 
are  plotted  in  Fig.  4.  The  oxygen  saturation  in  Fig. 
4(a)  follows  the  theoretical  Hill  curve28  reasonably 
well  with  a  mean  error  of  7.7%,  with  the  worst  accu¬ 
racy  close  to  zero  P02  and  the  higher  accuracy  when 
P02  is  above  80%  saturation.  For  P02  below 
20-mm  Hg,  oxygen  saturation  is  still  accurate  to 
within  15%,  with  a  low  standard  deviation  in  the 
images. 

With  variation  of  P02,  the  total  hemoglobin  con¬ 
centration  stayed  approximately  constant  [Fig.  4(b)] , 
with  a  mean  value  of  17.5  ±  2.1  jxM,  which  is  within 
97%  of  the  true  value  given  above,  and  water  exhib¬ 
ited  a  mean  value  of  94.2%  ±  8.3%.  Both  parameters 
show  some  cross  talk  at  P02  values  below  11-mm, 
which  is  possibly  unavoidable  owing  to  the  limited 
number  of  wavelengths  used  in  these  data.  Both  scat¬ 
ter  amplitude  and  scatter  power  stay  approximately 
constant  until  a  P02  of  1 1-mm  Hg,  beyond  which  both 
show  some  variation,  which  could  be  the  result  of 
cross  talk  between  the  two  parameters.  Above 
11-mm  Hg,  the  scatter  amplitude  has  value  of  0.92 
±  0.04  10-36(mm)6-1,  and  the  scatter  power  has  val¬ 
ues  of  1.49  ±  0.14.  The  reduced  scattering  coefficient, 
however,  stays  constant  throughout,  as  shown  for 
785  nm  in  Fig.  4(d),  with  a  mean  of  1.3 
±  0.03  mm-1. 


E.  Total  Hemoglobin  Validation 

The  final  experimental  homogeneous  data  set  involved 
varying  the  total  hemoglobin  while  keeping  oxygen 
saturation,  water,  and  scatter  parameters  constant. 
This  was  accomplished  by  use  of  a  similar  liquid  phan¬ 
tom,  with  1%  Intralipid  in  saline,  and  by  variation  of 
the  concentration  of  whole  blood.  The  hematocrit  level 
was  measured  by  a  clinical  co-oximeter,  yielding  1% 
blood,  which  is  equivalent  to  22-jxM  total  hemoglobin. 
The  blood  concentration  was  varied  from  0.2%  to  1% 
in  increments  of  0.2%,  and  the  amplitude  and  phase 
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Fig.  4.  Estimated  mean  values  are  shown  from  homogeneous  reconstructions  of  a  phantom  with  varying  oxygen  partial  pressure  (P02) 
of  the  solution  induced  by  addition  of  yeast.  The  oxygen  saturation  is  shown  in  (a),  along  with  the  theoretical  estimate  from  the  Hill  curve. 
The  total  hemoglobin  and  percent  water  are  shown  in  (b)  with  the  theoretically  estimated  values  of  18  mM  and  100%,  respectively,  and 
should  not  vary  with  changes  in  oxygenation.  The  scatter  power  and  amplitude  are  shown  in  (c)  and  should  not  vary.  The  reduced  scattering 
coefficient  at  785  nm  is  also  shown  in  (d). 


measurements  were  obtained  for  each  level.  After 
calibration,  the  spectrally  constrained  reconstruction 
was  applied  to  the  data,  and  the  recovered  mean  and 
standard  deviation  from  the  NIR  parameters  are 
plotted  in  Fig.  5.  The  theoretical  water  and  oxygen 
saturation  values  were  100%,  as  the  phantom  was  a 
liquid  oxygenated  solution. 

Figure  5(a)  shows  that  the  total  hemoglobin  fol¬ 
lowed  the  variation  in  blood  (%)  linearly  and  is  quan¬ 
titatively  accurate,  with  a  mean  percent  error  of 
6.2%.  Oxygen  saturation  stayed  constant  with 
change  in  blood  concentration  [Fig.  5(b)] ,  with  a  mean 
value  of  98.9  ±  0.6%.  The  same  trend  was  found  in 
water,  with  a  mean  value  of  98.2%  ±  1.5%.  The  scat¬ 
ter  amplitude  and  scatter  power,  shown  in  Fig.  5(b), 
are  also  independent  of  the  variation  in  blood  concen¬ 
tration,  with  a  scatter  amplitude  of  0.65  ±  0.01.  The 
scatter  power  had  a  mean  value  of  1.39  ±  0.08,  and 
this  agreed  well  with  the  estimated  1.4  derived  from 
van  Staveren  et  al. 19  The  reconstruction  converged  in 
4-6  iterations  for  the  different  concentrations,  and 
no  spatial  filtering  was  applied  in  the  reconstruction. 

4.  Discussion  and  Conclusions 

The  spectrally  constrained  direct  chromophore  and 
scattering  reconstruction  has  been  implemented  and 


validated  by  use  of  simulated  and  experimental 
frequency-domain  measurements.  The  results  from 
Subsection  3. A  show  improved  robustness  of  the  re¬ 
construction  to  increased  amounts  of  noise  in  the 
data.  The  frequency-domain  instrumentation  in  use 
typically  has  0.5%  noise  in  amplitude  and  0.5  deg  in 
phase,11  and  at  this  noise  level.  With  simulated  data 
from  a  homogeneous  phantom,  water  images  show  a 
reduction  in  standard  deviation  from  32%  to  10%  (as 
percent  of  the  mean)  in  going  from  the  conventional 
to  the  spectral  approach.  Even  at  a  5%  noise  level,  the 
spectral  approach  shows  recovery  of  the  parameters 
that  is  accurate  within  10%,  on  average,  with  a  sig¬ 
nificantly  reduced  standard  deviation  as  compared 
with  the  conventional  method.  This  insensitivity  to 
noise  is  due  to  the  use  of  multiwavelength  data  to¬ 
gether  with  the  spectral  constraints,  which  results  in 
a  reduction  in  the  number  of  unknowns,  making  the 
reconstruction  problem  less  ill  posed.  The  reduction 
in  noise  in  the  images  is  also  observed  in  the  homo¬ 
geneous  experimental  data  and  is  especially  signifi¬ 
cant  in  water  and  scatter  power.  The  NIR  parameter 
images  in  Fig.  2  show  the  suppression  of  the  artifacts 
in  the  images,  particularly  in  water  and  scatter 
power,  that  is  similar  to  the  trend  observed  in  the 
simulated  data. 
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Fig.  5.  Reconstructed  mean  and  standard  deviation  values  are 
shown  from  phantoms  with  varying  concentrations  of  blood  in 
which  the  total  hemoglobin  is  graphed  (a)  alongside  the  theoretical 
value  shown  by  a  dotted  line.  The  recovered  values  of  oxygen 
saturation  (%),  water  (%),  scattering  amplitude  (a)  and  scattering 
power  (b)  are  shown  in  (b),  all  of  which  are  not  expected  to  vary 
with  changes  in  total  blood  (%).  The  scatter  parameters  are  mul¬ 
tiplied  by  100  to  allow  them  to  be  displayed  on  the  same  graph  as 
the  oxygen  saturation  and  the  water. 


Water  is  an  important  measure  of  breast  physiol¬ 
ogy;  however,  its  quantitative  accuracy  from  NIR  to¬ 
mography  has  yet  to  be  validated.  In  the  past  several 
studies  used  a  fixed  water  content  in  tissue  (such  as 
assuming  30%-31%  fraction  in  tissue)  to  allow  for 
spectral  fitting  of  hemoglobin  levels29’30  or  used  ra¬ 
diological  data31  because  sufficient  wavelengths  were 
unavailable  and  the  cross  talk  between  water  and 
oxyhemoglobin  prevented  quantitatively  accurate  re¬ 
covery  of  the  water  absorption.  McBride  et  al.32  and 
Cerussi  et  al.33  have  shown  that  using  sparse  spectral 
image  data  from  a  subject  with  assumptions  about 
the  bulk  concentration  for  water  and  lipids  could 
have  up  to  15%  influence  on  HbT  and  oxygen  satu¬ 
ration  estimates.  In  recent  studies  water  has  been 
shown  to  have  significant  variation  with  breast  size 
in  normal  subject  studies,24  with  values  between  sub¬ 
jects  varying  significantly  from  10%  to  70%.  This 
large  range  was  observed  in  a  30-subject  population, 
as  shown  by  Cerussi  et  al.  33  and  values  of  21%  to  82% 
were  observed  in  a  26-subject  population  studied  by 
Srinivasan  et  al.24  These  large  numbers  suggest  that 
water  is  a  measure  of  the  extr avascular  space,  since 
the  vascular  space  is  clearly  less  than  2%  in  most 


breast  tissues.  Thus  water  yields  different  informa¬ 
tion  about  the  physiology  of  the  breast  than  does 
hemoglobin,  and  spatial  changes  are  expected  in  wa¬ 
ter  owing  to  different  content  in  the  fatty  and  glan¬ 
dular  tissues,  which  varies  with  the  composition  of 
the  breast.  The  change  in  water  content  during  the 
course  of  the  menstrual  cycle  has  been  followed  by 
Cubeddu  et  al. ,  showing  an  increase  in  water  in  the 
second  half  of  the  cycle34  for  one  patient.  Shah  et  al.35 
have  shown  an  increase  of  up  to  28.1%  in  water  in  the 
luteal  phase  in  a  single  volunteer,  and  Pogue  et  al.36 
have  shown  individual  variations  in  seven  subjects, 
with  a  mean  value  nearer  to  2.5%.  In  a  study  follow¬ 
ing  the  effect  of  neoadjuvant  chemotherapy37  in  a 
subject  with  a  palpable  adenocarcinoma,  using  opti¬ 
cal  white-light  spectroscopy,  water  showed  the  most 
dramatic  change,  dropping  67%  over  the  course  of 
three  treatment  cycles.  These  studies  have  shown 
important  trends  in  water  content.  However,  to  im¬ 
prove  the  clinical  utility  of  the  recovered  values,  it  is 
important  to  show  that  they  are  quantitatively  accu¬ 
rate  estimates  as  well.  The  results  shown  here  pro¬ 
vide  evidence  of  improvement  in  the  quantification  of 
water  by  use  of  the  spectrally  constrained  approach 
compared  with  the  conventional  method  of  fitting  for 
the  chromophore  and  scattering  parameters  from  op¬ 
tical  property  reconstructions.  Water  fraction  values 
obtained  by  the  new  approach  in  data  sets  that  follow 
variation  in  scattering,  oxygen  saturation,  and  total 
hemoglobin  agree  well  with  theoretical  predictions 
and  exhibit  reduced  noise  and  cross  talk  with  oxyhe¬ 
moglobin  compared  with  the  conventional  technique. 

Scattering  is  another  area  in  which  the  constraints 
from  Mie  theory  incorporated  into  the  reconstruction 
significantly  improve  the  quantification  of  the  scatter 
amplitude  and  power.  Recent  studies  have  used  scat¬ 
tering  to  study  structural  variations.  For  example, 
Poplack  et  al.38  showed  in  a  normal  cohort  of  23 
women  that  there  is  a  significant  decrease  in  the 
reduced  scattering  coefficient  (jju/)  at  785  nm  with 
increasing  body  mass  index  and  that  adipose  tissue 
was  less  scattering  than  glandular  tissue,  as  ex¬ 
pected.  The  same  trend  was  observed  by  Durduran  et 
al.29  in  a  subject  pool  of  52  volunteers.  Cerussi  et  al.33 
used  the  fit  of  p,s'  to  the  Mie  theory  approximation  to 
show  that  scatter  power  decreases  with  increasing 
body  mass  index  in  a  group  of  30  healthy  women. 
Pogue  et  al.36  showed  that  scattering  power  and  scat¬ 
ter  amplitude  could  successfully  separate  categories 
of  fatty  and  scattered  breasts  from  extremely  dense 
breasts  (p  <  10-4)  in  a  survey  of  39  women  with 
normal  mammography.  Since  the  risk  of  cancer  is 
strongly  correlated  to  the  breast  density,39  this  sep¬ 
aration  of  breast  densities  in  a  noninvasive  manner 
could  prove  very  useful.  Although  these  trends  may 
prove  to  be  promising  in  the  future,  their  quantitative 
accuracy  must  be  validated  in  phantom  studies. 
Hence  it  is  important  to  investigate  the  accuracy  and 
standard  deviation  in  these  parameters  to  fully  ex¬ 
ploit  the  NIR  information.  Results  in  Subsection  3. A 
show  the  improved  reconstruction  of  scatter  power 
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images  with  robustness  maintained  at  noise  levels  as 
high  as  5%  noise.  In  Fig.  2  the  recovered  mean  for  the 
scatter  amplitude  and  scatter  power  agree  well  with 
predicted  values  from  van  Steveran  et  al. 19  for  homo¬ 
geneous  experimental  data  (containing  1%  blood  and 
1%  Intralipid  in  saline).  A  change  in  scattering,  ob¬ 
tained  by  varying  the  Intralipid  concentration  in  a 
homogeneous  phantom  solution,  was  successfully  fol¬ 
lowed  (Fig.  3),  in  which  the  scatter  amplitude  varied 
linearly  and  the  scatter  power  showed  a  mean  of 
1.4  ±0.1  during  the  changes  in  concentration.  The 
cross  talk  between  scatter  parameters  and  chro- 
mophore  concentrations  in  minimal,  with  total  hemo¬ 
globin,  oxygen  saturation,  and  water  content  staying 
constant  through  this  change  [Figs.  3(c)  and  3(d)], 
except  at  the  highest  scattering  concentration  of  In¬ 
tralipid  (1.5%).  The  mean  of  p,s'  also  exhibits  a  linear 
increase,  and  p,s'  at  785  nm  doubles  as  the  Intralipid 
concentration  changes  from  0.75%  (|ju/  = 

0.89  mm1)  to  1.5%  (p,s'  =  1.7  mm-1),  shown  in  Fig. 
3(b). 

The  Hill  curve  relation  between  oxygen  saturation 
(S02)  and  P02  of  oxygen  was  followed  by  the  spectral 
method  in  the  graphs  of  Fig.  4,  which  showed  that 
S02  was  accurate,  with  a  mean  error  of  7.7%.  Total 
hemoglobin  estimates  remained  constant  through 
this  change  in  P02,  producing  a  mean  that  is  accurate 
to  97%  of  the  expected  value.  Water  content  also 
showed  this  trend,  with  a  mean  of  94.2%,  which  com¬ 
pared  well  with  the  predicted  value  of  100%.  Al¬ 
though  some  variation  was  found  in  the  scatter 
amplitude  and  scatter  power,  p,s'  at  785  nm  remained 
constant  with  change  in  oxygenation,  with  a  mean  of 
1.3  ±  0.03  mm-1.  The  standard  deviations  in  the  ox¬ 
ygen  saturation  images  remain  low,  even  at  lower 
oxygenation.  This  translates  into  a  successful  recov¬ 
ery  of  S02  in  malignancies  without  significant  noise 
or  cross  talk  between  scatter  parameters  and  deoxy¬ 
hemoglobin.  Total  hemoglobin  recovered  by  the  spec¬ 
tral  approach  was  separately  validated  by  use  of 
experimental  data  obtained  by  varying  blood  concen¬ 
tration  from  0.2%  to  1%.  Quantitatively  accurate  re¬ 
sults  with  a  mean  percent  error  of  6.2%  were 
obtained.  No  cross  talk  between  any  of  the  parame¬ 
ters  was  observed  during  this  variation,  as  shown  in 
Fig.  5(b)  in  which  oxygen  saturation,  water,  scatter 
amplitude,  and  scatter  power  remained  unchanged 
and  close  to  predicted  values. 

The  spectrally  constrained  approach  is  inherently 
robust  owing  to  the  addition  of  a  priori  spectral  be¬ 
havior.  It  requires  less  spatial  filtering,  whereas  the 
conventional  technique  benefits  from  a  mean  filter  to 
prevent  excessive  noise  in  the  images.  This  new  ap¬ 
proach  also  converges  faster  and  is  readily  extend¬ 
able  to  three-dimensional  models  as  well.  Use  of  data 
at  additional  wavelengths  can  easily  be  implemented 
without  much  computational  burden  in  the  inversion 
process.  A  partial  volume-type  reconstruction  may  be 
an  efficient  means  of  processing  a  large  number  of 
measurements.  There  is  certainly  some  cross  talk 
found  at  P02  values  of  oxygen  lower  than  11-mm  Hg 


or  under  extremely  high  scattering  conditions,  which 
can  probably  be  resolved  only  by  the  addition  of  more 
wavelengths.  Preliminary  simulations  have  shown 
this  to  be  true.  The  improved  quantification  of  and 
robustness  to  noise  of  the  reconstruction  shown  here 
for  homogeneous  measurements  is  currently  also  ap¬ 
parent  in  the  heterogeneous  data. 

As  the  use  of  NIR  tomography  expands,  spectrally 
constrained  reconstruction  should  add  considerable 
value  to  obtaining  quantitatively  accurate  estimates 
of  different  parameters,  particularly  water  and  scat¬ 
ter  power.  The  use  of  frequency-domain  measure¬ 
ments  allows  a  good  separation  of  chromophore  and 
scattering,  and  together  with  the  spectral  approach, 
we  obtain  reduced  cross  talk  between  the  parame¬ 
ters,  suppression  of  image  artifacts,  insensitivity  to 
noise  in  the  measurements,  and,  finally,  accurate 
quantification  of  the  NIR  parameters. 
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